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Abstract. The role of turbulence in various astrophysical settings is reviewed. 
Among the differences to laboratory and atmospheric turbulence we highlight the 
ubiquitous presence of magnetic fields that are generally produced and maintained 
by dynamo action. The extreme temperature and density contrasts and stratifications 
are emphasized in connection with turbulence in the interstellar medium and in stars 
with outer convection zones, respectively. In many cases turbulence plays an essential 
role in facilitating enhanced transport of mass, momentum, energy, and magnetic fields 
in terms of the corresponding coarse-grained mean fields. Those transport properties 
are usually strongly modified by anisotropics and often completely new effects emerge 
in such a description that have no correspondence in terms of the original (non coarse- 
grained) fields. (16 March 2011, Revision: 1.125). 



1. Introduction 

Astrophysical flows tend to be turbulent in the sense of being highly irregular. The 
study of astrophysical turbulence is important for several reasons. Firstly, turbulence 
needs to be taken into account when modeling most astrophysical systems. It can 
provide enhanced turbulent viscosity, turbulent heating, turbulent pressure, and leads 
to other effects, some of which can be non-diffusive in nature. Secondly, turbulence 
needs to be taken into account when interpreting observations of such systems. This 
is particularly evident in modeling line broadening and line asymmetries. Thirdly, 
astrophysical turbulence often spans an enormous range of length scales, allowing unique 
insights into the scaling properties of turbulence. 

In many text books various definitions of turbulence are suggested. However, none 
of them is quite without problems. Throughout this review, turbulence will remain a 
loosely defined property of flows that are highly irregular in space and time. 

Astrophysical turbulence as such is in principle not any different from ordinary 
turbulence. What is characteristic about it is the extremes in some parameters, e.g. 
huge Reynolds numbers, Prandtl numbers very different from unity, and, in some cases, 
strong density stratification and/or very high Mach numbers. Also, of course, the gas is 
often ionized and hence electrically conducting, so the interaction with magnetic fields 
cannot be neglected. As a rule, astrophysical flows tend to be magnetized spontaneously 
by self-excited dynamo action. 
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Contrary to laboratory and technical realizations of turbulence, where the driving 
often comes from the interaction with boundaries, astrophysical turbulence tends to be 
largely independent of explicit boundaries and is facilitated by intrinsic instabilities. 
Another difference between astrophysical and laboratory turbulence is the fact that, 
with very few exceptions, in situ observations are impossible and one has to rely on the 
radiative properties of the gas to infer velocity, temperature, and magnetic fields, for 
example. Yet another difference is that in some astrophysical flows the gas is extremely 
tenuous and close to collisionless, so the fluid approximation may actually break down. 
In some cases, multi-fluid descriptions are possible, for example when charged and 
neutral species move at different speeds, have different temperatures, or when positive 
and negative charge carriers, as well as dust particles need to be considered. However, 
quite often the multi-fluid description is then also not sufficient and it is better to employ 
more accurate techniques using, for example, particle-in-cell (PIC) methods or to solve 
the underlying Vlasov equations. This can be made more feasible by making use of 
the guide-field or gyro-kinetic approximations, where one averages out the azimuthal 
particle position around magnetic field lines. 

Astrophysical turbulence has been discussed in many excellent text books and 
reviews [H|2j[3jlH[5j|6l[71[8j|9j [TU] . In recent years, however, high resolution 
numerical simulations have become feasible and have added significantly to our 
understanding. Furthermore, the availability of three-dimensional codes has helped to 
make astrophysical turbulence a natural ingredient to many astrophysical investigations. 
The purpose of this review is to highlight recent progress in the field. We will focus on 
hydrodynamic and magnetohydrodynamic (MHD) aspects, but will try to keep the level 
of repetition with earlier reviews at a minimum. In particular, we shall not go in depth 
into aspects of dynamo theory, but refer instead to the recent review of Brandenburg 
& Subramanian [TTJ on recent progress and in particular on the nonlinear saturation of 
dynamos. 

2. Commonly used tools and conventions 

Throughout this review we assume some basic level of familiarity with commonly used 
tools and techniques in turbulence research. Here we only review the essentials in 
simplistic terms. 

2.1. Spectra 

Common tools include energy and helicity spectra, as well as structure functions and 
structure function exponents. These concepts become particularly useful if spatial 
homogeneity can be assumed. In simulations this usually means that one deals with 
triply periodic boundary conditions. Alternatively, one can apply these tools to just one 
or two periodic directions (for example convection in a domain with periodic boundary 
conditions in the horizontal directions). 
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In incompressible (or nearly incompressible) isotropic turbulence one usually defines 
the spectral energy per unit mass, 

E(k,t)= £ \u(k,t)\\ (1) 
fc_<|fc|<fc+ 

where k± = k ± 8k/ 2 mark a constant linear interval around wavenumber k, and the 
hat on u denotes the three-dimensional Fourier transformation in space. The spectral 
kinetic energy is normalized such that 

/ E(k)dk = l(u 2 ), (2) 
Jo 

where angular brackets denote volume averaging. This equation shows that the 
dimension of E(k,t) is cm 3 s" 2 , and E(k) can be interpreted as the kinetic energy 
per unit mass and wavenumber. 

In turbulent flows spectra remain in general time-dependent, so one is interested in 
spectra that are averaged over a sufficiently long time span. Such spectra can then also 
be compared with analytic predictions where statistical averaging is adopted instead. 

In strongly compressible flows one can also define the spectrum of kinetic energy 
per unit volume as 

E 2 (k,t)= £ \p^u\\ (3) 
fc_<|fc|<fc+ 

and the spectrum 

E 3 (k,t)= J2 lA^til 8 , (4) 
fc_<|fc|<fc+ 

which does not have a simple physical interpretation, except that Es(k,t) 3 ' 2 , integrated 
over k, has the dimension of an energy flux [12]. In strongly compressible (e.g. 
supersonic) flows these various spectra can become quite distinct. The closest agreement 
between spectra for subsonic and supersonic turbulence is achieved when using the 
quantity E 3 (k, t) [15] . 

In anisotropic turbulence it is useful to consider the spectral energy dependence 
along and perpendicular to the preferred direction of the turbulence, i.e. E(k±,k\\,t). 
Examples where this is important include rotating turbulence and turbulence in the 
presence of a strong magnetic field, but also inhomogeneous turbulence such as stratified 
turbulence and convection where one usually considers the spectral dependence on the 
horizontal wavenumber only. 

The kinetic helicity spectrum is defined as 

F(k) = ]T (w* -u + uj-u*), (5) 

fc_<|fc|<fc + 

where u) = V x u is the vorticity, and asterisks denote complex conjugation. The kinetic 
helicity spectrum is normalized such that 

/ F(k)dk = (w ■ u). (6) 
Jo 
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This kinetic helicity spectrum obeys the realizability condition, 

\F(k)\ < 2kE(k), (7) 

which is easily demonstrated by decomposing velocity and vorticity into positively and 
negatively polarized waves [H] H3]. Sometimes the helicity is defined with a 1/2 factor, 
just like the energy is. In that case the factor 2 in equation (|7|) would disappear. 

Equivalent concepts and definitions also apply to the magnetic field B, where 
one defines spectra of magnetic energy M(k), magnetic helicity H(k), and current 
helicity C(k), which are normalized such that / M(k) dk = (B 2 )/2fi , where /xq is 
the vacuum permeability, J H(k) dk = (A ■ B), and / C(k) dk = (J ■ B). Here, A is 
the magnetic vector potential with B = V x A and J = V x B/[1q is the current 
density. The magnetic helicity and its spectrum are gauge-invariant because of the 
assumed periodicity of the underlying domain. In that case the addition of a gradient 
term, VA, in A has no effect, because (VA • B) = (AV • B) = 0, where we have used 
the condition that B is solenoidal. Additional mathematical properties can be found 
in Ref. [15J. Magnetic helicity is an important quantity, because it is conserved in the 
limit of vanishing magnetic resistivity and in the absence of boundary losses. Another 
similarly conserved quantity is the cross-helicity, (u ■ B). Its sign indicates whether 
Alfven waves travel preferentially parallel or antiparallel to the local magnetic field. 

2.2. Turbulent cascade 

The energy-carrying scale is often defined as the scale if = 2ir/kf, where kf is the 
wavenumber where the energy spectrum peaks. It is close to the integral scale l\ = 2ir/ki, 
where fcf 1 = / k^E^k) dk/ J E(k) dk. 

Turbulence is driven either by some explicit stirring or by some type of instability. 
Explicit stirring is frequently used in direct numerical simulations (DNS) and large eddy 
simulations (LES). Here, DNS means that one considers the original equations with the 
proper diffusion term, as opposed to other schemes such as LES that are motivated 
by numerical considerations and limited resolution. An astrophysical example is the 
driving accomplished by supernova explosions in the interstellar medium within each 
galaxy. Examples of instability-driven turbulence include Rayleigh-Benard convection, 
the magneto-rotational instability (MRI), and shear flow instabilities with inflection 
points resulting from rigid surfaces such as the accretion disc near the surface of a 
neutron star. In the latter case the domain is obviously no longer periodic. 

The driving usually occurs over a certain range of length scales around the 
wavenumber kf. The nonlinearity of the hydrodynamic equations produces power on 
progressively smaller scales (larger wavenumbers) . Qualitatively, this leads to a cascade 
of energy from large to small scales until energy is dissipated at scales corresponding to 
the wavenumber k^. The range of wavenumbers between kf and k^ is called the inertial 
range. An important quantitative property of turbulence is the approximate constancy 
of spectral energy flux e throughout the inertial range, where e has dimensions cm 2 s -3 . 
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Making the ansatz 



E(k) = C K e a k\ 



(8) 



where Ck is the Kolmogorov constant, the values of the exponents a and b are determined 
by matching the dimensions for length (cm) and time (s) as follows: 3 = 2a — b and 
2 = 3a, respectively. This yields a = 2/3 and b = -5/3, so E(k) = C K e 2/3 k~ 5 / 3 . 

The length of the inertial range can be calculated by assuming that E(k) is finite 
only in the range kf < k < k&. Thus, u rms and e are given by the two integrals 



which are just the normalization condition of E(k) and the definition of the energy 
dissipation, respectively. Here, v is the kinematic viscosity. Eliminating e, and writing 
the result in terms of the Reynolds number yields 



Thus, the length of the inertial range scales with the Reynolds number like k^/kf f» 
Re 3/4 ; see, e.g., Ref. [16]. 



In astrophysics one often deals with extremely large Reynolds numbers, and hence 
an extremely broad inertial range. This is in practice not possible to simulate. However, 
many aspects of interest are independent of the length of the inertial range. Those that 
are not exhibit a well-understood scaling with Reynolds number. This is the main reason 
why it is at all possible to attempt simulating astrophysical systems on the computer! 

2.3. Taylor hypothesis and one- dimensional spectra 

In laboratory and atmospheric turbulence, for example, one usually measures time series 
which allow only one-dimensional spectra to be determined. This involves making the 
Taylor hypothesis, i.e. the assumption that the temporal power spectrum, u(uj), can be 
associated with a spatial one, u(k), via u = Uok. Here, Uq is the mean flow at the 
location of the detector. 

It is important to realize that one-dimensional spectra can differ from the fully 
three-dimensional spectra that are normally considered in numerical simulations of 
turbulent flows. The two agree only in regions of the spectrum where one has power 
law scaling, i.e. where E(k) ~ k n with some exponent n, This is evidently not the case 
near the dissipation subrange and near the sub-inertial range at small wavenumbers. 
This is probably the main reason why spectra from high resolution DNS show a 
significantly shallower spectrum just before the dissipative subrange than the one- 
dimensional spectra obtained using the Taylor hypothesis, where a shallower part in 
the spectrum is essentially absent. In the following we briefly explain this difference 




(9) 




(10) 
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Consider the case of a one-dimensional spectrum obtained by Fourier transforma- 
tion over the z direction. To relate this to the three-dimensional spectrum, we average 
over the remaining x and y directions. Thus, we compute for k z > 

E 1D (k z )= II \u(x,y,k z )\ 2 ^. (12) 

J J ±J X ±Jy 

Next, using Parseval's relation for converting the averaging in real space to an integration 
in spectral space, we can write 

where we have assumed that |ilt| 2 is statistically axisymmetric, i.e. independent of the 
azimuthal angle about the k z axis. Next, we use k 2 = k 2 — k 2 z to replace the k r dk r 
integration by one over kdk in the range from k z < k < oo, i.e. 

E 1D {k z ) = 2tt [°°\u(k)\ 2 kdk = [°°^-dk, (14) 

J k z J k z K 

where we have used the fact that the three-dimensional spectrum can also be written as 
E(k) = 2nk 2 \u(k)\ 2 , where we have assumed averaging over full shells in wavenumber 
space. Thus, we see that one-dimensional spectra, En)(k), are related to the fully 
three-dimensional spectra, E(k), via integration, or via differentiation for the reverse 
operation, i.e. 

E m (k) = r^p-dk' and E(k) = -k^^. (15) 
Jk k' dk 

We reiterate that, if one of the two spectra were a pure power law, the other one would 

also be a pure power law. However, this assumption breaks down near kf and k&. We 

mention this aspect here, because one of the unexpected results obtained from a number 

of simulations over the last decade is a strong departure from the Kolmogorov k~ 5 ^ 3 slope 

near k^, where the spectrum can be substantially shallower [TRl \T7\ ITT?] . This is now 

known as the bottleneck effect [20J and was first noticed in atmospheric turbulence {21] . 

It is by far not as marked in one-dimensional spectra as in three-dimensional spectra 

from recent high-resolution DNS |17j . 



E 1D (k z ) = J J \u(k x , k y , k z )\ 2 dk x dk y = 2vr J 



2.4- Intermittency 

The scaling of velocity differences over fixed distances is different in different locations. 
The flow is therefore said to be intermittent. A related property is that the scaling of 
the structure functions, 

S p (r,t) = (\u(x + r,t)- u(x,t)\ p ) , (16) 

with distance r = \r\ deviates from the scaling r p / 3 for all moments p ^ 3, for both 
parallel (r parallel to u) and transverse (r perpendicular to u) structure functions. This 
property is quantified by the structure function exponents, ( p , which denote the slopes 
in graphs of In S p (r,t) with lnr. The averaging, denoted by angular brackets, is here 
taken to be over the full volume. 
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In practice, approximate scaling can only be identified in a rather limited range 
of lnr. Analytic theory predicts £3 = 1; see, e.g. Ref. [61 |8]. This property is often 
used to improve the accuracy in the determination of ( p for p 7^ 3 from numerical or 
experimental data by plotting In S p (r, t) versus In Ss(r, t). This procedure is referred to 
as extended self-similarity or ESS [22J. 

Intermittency is linked to the property that the ( p deviate from a linear 
dependence on p. A completely non-intermittent behavior would mean ( p = p/3. A 
phenomenological relation that describes the behavior observed in experiments and 
simulations is given by the She-Leveque relation [23J 

P/31 



4 = | + c 



(17) 



where C is interpreted as the co-dimension of the dissipative structures. For weakly 
compressible or incompressible turbulence the dissipative structures are one-dimensional 
tube-like structures, so the co-dimension is C = 2. Under compressible conditions the 
dissipative structures tend to become two-dimensional sheet-like structures, so C = 1 
[24J, which is also borne out by simulations of highly supersonic turbulence (251 EE]. 
Sheet-like dissipative structures are also expected in hydromagnetic turbulence, where 
these structures correspond to current sheets. In that case one expects the same scaling 
as for supersonic turbulence [27]. However, in incompressible hydromagnetic turbulence 
with constant density p = p , the relevant structure functions are based on the so-called 
Elsasser variables z — u ± B/ y/JIoPo- In that case, analytic theory predicts that the 
mixed third-order longitudinal structure functions of Politano & Pouquet [28] , 

S%(r) = (5zf(r)[6z±(r)] 2 ), (18) 

scale linearly with r = \r\. Here, 8z ± (r) = z ± (x + r) — z ± (x), 5zu~(r) = 5z ± ■ f, and 
f = r/r is the unit vector of r. 

Simulations tend to give slightly different scalings for the longitudinal and 
transverse structure functions. This may be a consequence of different cascade speeds for 
longitudinal and transverse velocity increments [29], but it may also just be an artifact 
of insufficient resolution and may go away at larger resolution, as indicated by recent 
simulations at high numerical resolution [T3] . 

The assumption of the constancy of the spectral flux is well confirmed, but the 
correlation between energy injection and energy dissipation displays significant scatter. 
This is mostly because the spectral flux fluctuates significantly in time and there is 
some delay before the spectral energy has reached the dissipation scale. By taking into 
account the appropriate delay the scatter can be significantly reduced [30]. The energy 
flux at large scales is characterized by 

e = C e ul D /L, (19) 

where C e ~ 0.5. It is customary to define the length scale as L = 37r/4fcf, so in terms of 
kf and u^ ms = 3u^ D we can then write 

0.04 k f u 3 ims . (20) 
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This formula will be useful later in connection with turbulence in interstellar and 
intergalactic media. 



3. Sites of astrophysical turbulence 

The following discussion is concerned mainly with observations and simulations covering 
a range of astrophysical settings where turbulence occurs. In some cases strong 
theoretical evidence is used to argue for the existence of turbulence, as for example 
in accretion discs where turbulence has not yet been observed explicitly [31 J. 



3.1. Solar wind 

The gas above the visible surface of the Sun is not in hydrostatic equilibrium. Instead, 
because of geometrical constraints and because of a gravitational potential inversely 
proportional to the radial distance, there is the possibility of a critical point, where the 
radial velocity equals the sounds speed. The theory of such flows was first understood 
by Parker [32] in 1967 and is now explained in a number of text books on compressible 
flows or on astrophysical fluid dynamics [331 El]- Other transonic flows of this type 
include those through a Laval nozzle, as well as Roche-lobe overflow between binaries, 
astrophysical jets from accretion discs, and Bondi accretion. In the case of the Sun the 
gas reaches speeds of around 400 km s _1 in the equatorial plane and 800 km s _1 at higher 
latitudes [35]. The solar wind is turbulent and fluctuates between 300 and 800 km s -1 
on time scales ranging from seconds to hundreds of hours [7J. 

In the case of the solar wind, spectral information can be obtained under the Taylor 



hypothesis that was discussed in §2.3, Using this hypothesis the following properties 
have been inferred: 

• An approximate k~ 5 ^ 3 energy spectrum both for velocity and magnetic field [36J. 

• Below the ion Larmor radius a steeper spectrum (between k~ 2 and k~ A ) is found for 
the magnetic field [37]. In view of theoretical expectations the transition to a k~ 7 l 3 
spectrum for the magnetic field together with a k~ x ^ 3 spectrum for the electric field 
is particularly interesting [3HI |3J5] ; see Figure [TJ 

• Finite magnetic helicity (negative in the northern hemisphere and positive in the 
southern hemisphere), possibly with a k~ 7 ^ 3 spectrum |40j . 

• Finite cross helicity of positive sign, indicating outward travelling waves [41]. 

• Decay of turbulence with distance and evidence for additional heating [7J [3J)J H21 H3J . 

A possible connection between a k~ 7 ^ 3 tail in the energy spectrum at small scales 
(below the scale of the ion Larmor radius) and so-called electron MHD as a model 
for collisionless plasmas such as the solar corona and the Earth's magnetosphere has 
been discussed [44J. A similar slope has now also been seen in simulations using the 
gyrokinetic equations [38J. These equations emerge from the Vlasov equations for a 
collisionless plasma by averaging over the azimuthal angle of the gyrokinetic motions 
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Figure 1. Spectra of electric and magnetic fields from a gyrokinetic simulation [38] 
(left) compared with those obtained from the Cluster spacecraft [39] (right). Note the 
approximate fc" 5 / 3 spectrum below the Doppler-shifted inverse proton Larmor radius 
and an approximate fc -7 / 3 spectrum for the magnetic field (solid/blue on the left and 
light shade/green on the right) between the Doppler-shifted inverse proton and electron 
Larmor radii (in the right hand plot referred to as f pp and f pe , respectively), followed 
by a steeper dissipation subrange. Above the inverse Doppler-shifted electron Larmor 
radius the electric field spectrum develops a shallower subrange consistent with A; -1 / 3 
(dashed/red on the left and black on the right). Courtesy of Gregory Howes [38] . 
as well as Fouad Sahraoui and Melvyn Goldstein [3U], and Copyright (2010) by The 
American Physical Society. 

[35] . Let us also mention here the possibility of obtaining spectra steeper than /c~ 7 / 3 
using electron MHD when equipartition between kinetic and magnetic energies is not 
satisfied [16], or when compressible effects are included [17] . 

In view of our discussion in § |2.3[ it should be noted that near the break point 
where the spectral index changes, the spectra inferred by using the Taylor hypothesis are 
not exactly representative of the three-dimensional spectra obtained from simulations. 
However, in view of other general uncertainties, the changes in the spectral slopes are 
probably sufficiently weak to be ignorable. 

3.2. Solar convection 

The visible surface of the Sun is the photosphere, from where photons can reach the 
Earth in a direct path. Deeper inside the Sun the gas is opaque and photons are 
continuously absorbed and re-emitted, following approximately a diffusion-like process. 
At the surface, the Sun exhibits a granular pattern that can already be seen with 
small amateur telescopes. The pattern is irregular and changes on a time scale of 




Figure 2. Comparison between a granulation pattern from a simulation with 12 
km grid size (left), an observed granulation pattern from the Swedish 1-meter Solar 
Telescope at disk center (middle), and the simulated one after convolving with the 
theoretical point spread function of a 1 meter telescope. The simulation images are 
for wavelength integrated light intensity while the observed image is for a wavelength 
band in the near UV. The image was taken on 23 May 2010 at 12:42 GMT with image 
restoration by use of the multi-frame blind de-convolution technique with multiple 
objects and phase diversity 48 . Courtesy of V. M. J. Hcnriques and G. B. Scharmer. 



around 5 minutes. The horizontal pattern size is 1-2 Mm. Here and elsewhere we use 
1 Mm = 1000 km as a convenient length scale. The visible granulation is just a thin 
layer on top of a 200 Mm deep convection zone. The convection zone covers the outer 
30% of the Sun by radius. The inner 70% are convectively stable. This region is referred 
to as the radiative interior. 

Its overall dynamics can be understood through simulations and turbulence theory 
(i.e. mixing length theory) [S]. Excellent agreement between observed and simulated 
granulation patterns has been obtained; see Figure [2| By calculating diagnostic spectra 
in the visible light and comparing with observations one can determine the abundance 
of chemical elements \$$\ \5U\ lol] . The chemical element abundances are important for 
determining the opacity of the gas which, in turn, determines the radial structure of the 



Sun. This will be discussed in more detail in § 6.10 

From the viewpoint of turbulence theory, this type of convection is special - not so 
much because the Rayleigh number is extremely large (~ 10 30 ), but mainly because the 
density and temperature stratifications are extreme, covering 6 orders of magnitude 
of change in density and a factor of 300 in temperature. This huge stratification 
implies that the turbulence characteristics become strongly depth-dependent. It has 
long been anticipated that the energy-carrying scale varies with depth in such a way 
that it is proportional to the local pressure scale height, H p . The pressure scale height 
is proportional to the temperature and varies from about 200 km at the top of the 
convection zone to about 60 Mm at the bottom. The typical correlation time of the 
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x [Mm] 

Figure 3. Vertical velocity near the surface from a solar convection simulation. The 
right hand panels show part of the domain at two different times. 

turbulence is expected to be proportional to the local turnover time, H p /u Tms , where 
u rms is the rms velocity of the turbulence. Estimating the convective energy flux as 
-Fconv ~ P u rmsy we expect u Tms to vary by a factor of 100 from about 4km/s at the top 
of the convection zone to about 40m/s at the bottom. Thus, the turnover times vary 
by more than 4 orders of magnitude, from minutes at the top of the convection zone to 
about a month at the bottom. 

A general difficulty in carrying out simulations of the deep solar convection zone 
is the long Kelvin-Helmholtz time in deeper layers. The Kelvin-Helmholtz time can 
be defined as the ratio of thermal energy density to the divergence of the energy flux 
or (operationally more convenient) as the total thermal energy above a certain layer 
divided by the solar luminosity. This time scale determines the thermal adjustment 
time and can be rather long. However, by preparing initial conditions such that the 
mean stratification as well as the fluctuations are close to those in the final state, the 
difficulty with long time adjustment times can be alleviated. 

Figure [3] shows an example from radiation hydrodynamics simulations of the 
horizontal pattern of the vertical velocity near the surface, and Figure [4] the same at a 
depth of about 4 Mm. One sees clearly that the number of cells has decreased and that 
the horizontal scale of the cells changes from about 2 Mm near the top to about 10 Mm 
at a depth of about 3 Mm. This illustrates two important properties: (i) The horizontal 
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Figure 4. Same as Figure [3| but for a deeper level, about 4 Mm away from the 
surface. 

cell size below the surface is typically a few times the distance from the surface (which 
really reflects that it is several times the local pressure and density scale heights), and (ii) 
the structure size increases so rapidly with depth that even using the concept of " cells" 
may be misleading. One also sees that the typical cell life time changes rapidly with 
depth. Near the surface the cell pattern shows some clear changes after only 3 min, while 
at 4 Mm depth the changes remain more limited even after 20 min. The scales of the 
patterns and their rates of change are thus generally consistent, at a semi-quantitative 
level, with assumptions that have generally been made in simplified (mixing-length) 
models of convection: 

• The energy- carrying scales of the turbulence are of the order of the local vertical 
pressure scale height, H p = |Vlnj9| _1 . 

• The turbulence varies on time scales comparable to the turnover time defined as 

Hp/ M rms . 

One should, however, not conclude that the numerical results 'confirm' a scaling with 
the pressure scale height. Mass conservation really involves the density scale height 
rather than the pressure scale height, and the main reason that analytical theories of 
convection have generally tended to avoid using the density scale height is that, because 
of a rapid change of the degree of hydrogen ionization there is a narrow layer close to 
the surface of stars where the density scale height may tend to infinity. 
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Many of the qualitative expectations from mixing length theory are borne out by 
simulations. This also includes the scaling of velocity and the temperature fluctuations 
with convective flux and hence with depth. Indeed, one finds that the convective energy 
flux (or enthalpy flux), F conv , is proportional to the negative specific entropy gradient. 
Velocity and temperature fluctuations scale like F^ v and F^ v , respectively; see Fig. 11 
in Ref. [52]. 

Early ideas about distinctively different modes of convection at different scales 
are mostly due to differences in observational techniques rather than real physical 
differences in the convection. Supergranulation, for example, refers to a convection 
pattern with a horizontal scale of about 30 Mm, which is seen in Dopplergrams measuring 
the line of sight velocity. When plotting the horizontal velocity amplitude as a function 
of horizontal size the supergranulation scales appear to be just a part of a rather 
featureless power law extending over many orders of magnitude in size [53]. Banana 
cells, on the other hand, refer to a theoretically expected pattern of convection in deeper 
layers. This expectation is based on the Taylor-Proudman theorem [54J, rather than an 
observationally established fact, but it remains a pronounced feature of convection in 
rotating shells between ±30° latitude [551 EE] ■ 

3. 3. Other effects of solar turbulence 

There are a number of properties that occur on scales that are larger than the energy- 
carrying scale. These properties include: 

• The angular velocity varies by about 30% in latitude (slow at the poles and fast at 
the equator) with approximate solid body rotation below the convection zone and 
a general deceleration in the outer 5% of the solar radius [57] . 

• There is a large-scale magnetic field exhibiting a 22 year cycle (11 years for the 
sunspot number) and a statistical antisymmetry of the radial field with respect to 
the equator (Figure [5]). 

• The solar surface exhibits a magnetic field that is strongest inside sunspots, where 
it is seen through Zeeman splitting. 

• Magnetic and current helicity with strong fluctuations, but well-defined averages: 
negative in the northern hemisphere and positive in the southern hemisphere; see, 
e.g., Fig. 1 in Ref. [5S]. 

In addition to the convective motions of the Sun, there are coherent wave patterns 
that correspond to discrete frequencies in wavenumber and frequency space. Using a 
technique called helioseismology [59] EQl EU [62], the information contained in these 
modes can be used to infer the depth dependence of sound speed and hence the radial 
dependence of the temperature of the Sun. 

Helioseismic constraints of the core temperature were important in pinning down 
the origin of the low observed neutrino flux from the Sun [H31 EH ES] in terms of neutrino 
oscillations, i.e. the Mikheyev-Smirnov-Wolfenstein effect [66| IS"?]. 
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Figure 5. Longitudinally averaged radial component of the observed solar magnetic 
field as a function of cos(colatitude) and time. Dark (blue) shades denote negative 
values and light (yellow) shades denote positive values. Note the sign changes both in 
time and across the equator. Courtesy of R. Knaack. 



Solar rotation lifts the degeneracy of modes with different azimuthal order and 
allows a determination of the dependence of the internal angular velocity on radius and 
latitude [57]. Rotation also causes the convection pattern to propagate in a prograde 
direction [68] . 

At the equator, the Sun rotates with a period of about £ rot = 26 d, but at the poles 
it spins about 30% slower. This is referred to as differential rotation. The angular 
velocity is Q = 2n/t TOt , but in helioseismology one often talks about the rotation rate, 
Q/2tt, which is measured in nHz. The equatorial value at the surface is 452 nHz. The 
radiative interior is found to rotate rigidly [57] . The interface between the differentially 
rotating convection zone and the rigidly rotating radiative interior is referred to as the 
tachocline [691 H3- 

3.4- Interstellar turbulence 

The gas between the stars can be observed in absorption or emission both at infrared and 
radio wavelengths. The line of sight velocity component can be determined by Doppler 
shifts of spectral lines; see, e.g., Ref. [7_T]. There is a general power law scaling of 
velocity amplitudes and velocity differences with geometrical scale j7_U [721 E3] . Velocity 
dispersions scale with size to a power of about 0.4 from sub-parsec scales to scales of the 
order of about 1 kpc; see Fig. 1 of Ref. [72]. The velocity scaling is practically the same 
in regions with varying intensity of star formation, indicating that the velocity scaling is 
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inertial, and driven mostly by energy input at large scales, rather than a result of direct, 
local driving by on-going star formation [TH [75j [76] . Direct evidence of turbulence on 
small length scales (~ 10 12 cm) in the ISM comes from radio scintillation measurements 

[771 [78]. 

Galaxies such as our own have typical radii of R ~ 15 kiloparsecs (kpc). Here, 
1 kpc = 3 x 10 21 cm is used as a convenient length scale. The density decreases rapidly 
away from the midplane with a typical density scale height of H « 70 pc. Near the 
midplane of a typical galaxy the 3D rms turbulent velocities are around 15kms _1 . This 
implies a typical turnover time, H p /u TTns , of around 5Myr (megayears). 

An important aspect is the occurrence of supernovae, which mark the death of 
massive stars and provide a significant energy release into the interstellar medium 
through thermal energy and momentum injection. Traces of supernovae are seen as 
supernova remnants, which give a qualitative idea about the nature of interstellar 
turbulence. 

Supernova explosions contribute about -Esn = 10 51 erg per explosion. With about 
20 supernovae per million years per kpc 2 estimated for the solar neighborhood this 
corresponds to an energy injection per unit area of 

J e SN dz « 2010 51 erg/(310 13 s x 9 10 43 cm 2 ) « 7 lO^ergcnT 2 s -1 . (21) 

This is almost two orders of magnitude more that what is required to sustain the 
turbulent energy dissipation per unit area and time, which, from equation (19), may be 
estimated to be 

J edz « 0.5pw 3 D « 10~ 24 gcnT 3 (lO^ms^ 1 ) 3 « 10" 6 erg cm" 2 s _1 , (22) 

where the mean density of the interstellar medium is p ~ 2 10 _24 gcm -3 and the one- 
dimensional rms velocity is Uxd ~ lOkm/s = 10 6 cms _1 . This is also in good agreement 
with simulations [79]. A visualization of density and magnetic field strength in such a 
simulation is shown in Figure [6] 

The linear polarization properties of synchrotron radiation can be used to infer the 
magnetic field both along the line of sight via Faraday rotation and perpendicular to it 
through the polarization plane projected onto the sky [801 El E2]. The field strength 
is typically around 5 /zG in the solar neighborhood of our Galaxy, but it can be several 
mG in the galactic center [83| [M] • For many spiral galaxies large-scale magnetic fields 
have been found. In many of them the magnetic field is approximately axisymmetric 
and symmetric about the midplane [85] . 



3.5. Accretion discs 

Accretion discs are disc-like structures through which gas gradually spirals toward a 
central massive object while converting potential energy into kinetic and magnetic 
energies that are dissipated and radiated away. This conversion is believed to be of 
turbulent nature and may be driven by the magneto-rotational instability [311 [86] . An 
alternative mechanism for disk dissipation is that the disk functions as a self-regulating 
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Figure 6. Two-dimensional slices through a three-dimensional simulation domain 
of supernova-driven turbulence in the interstellar medium showing the vertical 
distribution of the density (left) and magnetic field (right). Note the appearance of 
supernova remnants in density and magnetic fields as well as an overall concentration 
around the midplane at z — 0. Courtesy of Miguel de Avillcz 79J. 



buffer. As long as the disk accretion towards the central object is smaller than the 
rate of mass in-fall onto the disk from the surrounding nebula, the mass density of the 
disk increases. When the surface density reaches a level sufficient for gravitationally 
driven instabilities to develop, spiral waves starts to grow, develop into spiral shocks, 
and dissipation in the shocks then enhances the disk accretion enough to balance the 
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rate of in-fall onto the disk [87] . 

In order to allow material to spiral inward at a mass accretion rate M, half of 
the orbital potential energy must be converted viscously and resistively into heat and 
radiation. Therefore the total (bolometric) luminosity of an accretion disc is [88] 

GMM . . 

L = -Yn~> (23) 

where M is the mass of the central object and i?; n is the inner radius of the accretion disc. 
Obviously, the further the disc stretches toward the central object, i.e. the smaller the 
value of i?i n , the more efficient the energy conversion will be. Discs around black holes are 
most efficient in this respect, because here the innermost stable orbit is 1-3 Schwarzschild 
radii, i.e. (2-6) xGM/c 2 , where c is the speed of light. Thus, L = 0.1 x Mc 2 , which 
constitutes a much more efficient conversion than nuclear fusion, where the efficiency is 
only 0.007 x Mc 2 . Here we have used for M the rate of hydrogen burning [88]. Note 
that the factor 0.007 comes from the relative mass difference between a helium atom 
(4.0026) and four hydrogen atoms (1.0078). 



3.6. Turbulence in galaxy clusters 

Galaxies themselves tend to cluster on Mpc scales. There are typically around 10 4 
galaxies in a cluster, but some clusters can be substantially smaller. All clusters are 
generally strong X-ray emitters, but some are also strong radio-emitters resulting from 
synchrotron emission in the presence of magnetic fields. 

Typical temperatures are around 10 8 K corresponding to a sound speed of around 
1000 km s -1 . The implied velocity dispersion is also of that order, as expected when the 
system is in approximate Virial equilibrium. With typical length scales on the order 
of the density scale height, H p = 100 kpc, the turnover time is 100kpc/(1000km/s) = 
0.1 Gyr. This would also be the typical decay time of the turbulence in the absence of 
mechanisms driving the turbulence. 

Mechanisms for driving such turbulence include mutual encounters of clusters 
[89] 190] . Given that only a fraction of all galaxy clusters also have strong radio halos 
[DTj , one may speculate that these clusters have undergone a recent encounter or merger 
with another cluster within the last few gigayears. Obviously, in this scenario one would 
just have decaying turbulence between encounters. In the context of galaxy clusters this 
subject has been studied by various groups [52J E21 US]- Another mechanism that has 
been discussed in the literature is the driving by individual galaxies moving through the 
cluster and producing a turbulent wake behind them [H51 ESI I9"T] . 



3.7. Decaying turbulence in the early universe 

Various mechanisms for the generation of "primordial" fields have been proposed [98J. 
One problem is that the predicted magnetic field strengths are extremely uncertain. 
Another general problem is the small length scale of such fields. For example, after 
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Wave number 

Figure 7. Magnetic energy spectra at different times (increasing roughly by a factor 
of 2). The curve with the right-most location of the peak corresponds to the initial 
time, while the other curves refer to later times (increasing from right to left). Note the 
temporal growth of spectral magnetic energy at wavenumbers to the left of the peak 
and the associated propagation of spectral energy to successively smaller wavenumbers, 
i.e. to successively larger scales. Adapted from Refs. [1011 1102] . 

the electroweak phase transition, about 10 -10 s after the Big Bang, the horizon scale 
was around 3 cm. Magnetic fields generated during such a phase transition may possess 
magnetic helicity, but this is also rather uncertain [99] . However, during the subsequent 
decay of a helical field, energy is transformed to larger scale by an inverse cascade of 
magnetic helicity \100\ 1101] . Figure [7] shows the evolution of the resulting magnetic 
power spectrum at different times from a direct numerical simulation of the relevant 
hydromagnetic equations [102J. Simulations have demonstrated that turbulence decays 
in power law fashion with the total energy being proportional to t~ n , where n = 0.5 
for maximally helical fields and n = 1 for non-helical fields j!02j . By comparison, 
nonhelical fluid turbulence leads to n = 1.2 |103[ 1104] . As argued by Biskamp & Miiller 
|105] . helical fields may be more typical than non-helical ones. Of course, the magnetic 
field generated in rotating bodies (stars and galaxies, although neither is relevant to the 
early Universe) tend to be helical, but of opposite sign in the two hemispheres, so the net 
magnetic helicity would cancel to zero. On the other hand, the helical contribution of a 
field generated at an early phase transition will decay more slowly than the non-helical 
contribution, and so the relative importance of the helical fields will grow with time. 

The question of the decay law of helical MHD turbulence is still not fully settled. 
It is generally believed that the magnetic energy, Em, follows a power law decay, i.e. 
Em ~ t~ n , but proposals for the value of n range from 2/3 to 1/2, depending essentially 
on the assumptions made about the evolution of the typical length scale L of the energy- 
carrying motions. If one assumes L to be controlled by a resistive evolution of magnetic 
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helicity, Hm, i-e. 

dH M 277 



Am, (24) 



dt L 2 ' 

then, for a power law evolution of Hm we have L ~ t 1//2 , and with £"m = Hm/L and 
i?M ~ const we find [106J, 

E M ~ t- 1 ' 2 . (25) 

On the other hand, if one discards resistive effects, and assumes instead that the decay 
is controlled by inertial range turbulence, i.e. 

r3 ^3/2 ^5/2 

(26) 



dE M _ U 3 E M E M 



dt L L H M ' 

then, after integration, we obtain [105J 

E M ~ t" 2/3 , (27) 

together with L ~ t 2 / 3 . Note than in either of the two proposals one has assumed that 
Hm ~ LEm ~ const. However, in the former approach Hm is not assumed to be constant 
exactly, but to decay resistively like Hm ~ t~ 2r]t l L2 , which implies a corresponding speed- 
up of the decay of Em and hence an increase of n from 1/2 to 1/2 + 2-qt/L 2 . This may 
explain why simulations at finite rj |105l 1107] suggest exponents close to n = 2/3. This 
question needs to be followed up again in future at higher resolution, but simulations 
at moderate resolution have confirmed the idea of a correction factor proportional to 

t -2r,t/L* in the decay of Em fXQ6] , 

It is still unclear whether such primordial magnetic fields would have a detectable 
effect on the polarization signal of the cosmic background radiation and whether 
significant fields may have been present when the first stars or galaxies were formed. 
These questions are subject to current investigations |108l 1109] . Another subject under 
active investigation concerns the production of gravitational waves from the Maxwell 
stress associated with primordial magnetic fields |110[ II 1 11 1112[ 1113] . 



4. Theoretical studies of turbulence 

4-1. Incompressible turbulence 

Most turbulence research is restricted to incompressible turbulence, in which case the 
Navier-Stokes equations take the form, 

^ = - Vp + f + z/VV V-u = 0. (28) 

Here, D/Dt = d/dt + u ■ V denotes the advective derivative. It is this term that 
constitutes the important nonlinearity of the Navier-Stokes equations. In order to 
understand the nature of the nonlinearity it is useful to make use of the vector identity 
u ■ Vm = uj x u + ^Vw 2 , with uj = V x u. Thus, we have 

^=ttxw-Vp + / + uV 2 u. (29) 
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Owing to incompressibility, we have p = const, and only the reduced pressure, 
V = p/p+^ u2 ^ enters in equation (29). However, because of the solenoidality constraint, 



V • u = 0, the pressure gradient also constitutes a quadratic nonlinearity of the form 
p = V 2 V ■ (ux uj + f). (30) 



This relation follows directly from equation (29) after taking its divergence and noting 
that V • du/dt = V • V 2 u = 0. 

4-2. Compressible fluid dynamics 

In the compressible case, the Navier-Stokes equation can be written in the form 

p^ = -Vp + F + V-r, (31) 

where r = 2puS is the stress tensor, here assumed to be proportional to the kinematic 
viscosity v and the traceless rate of strain tensor, S, whose components are 

Sij = \ {uij + Uj t i) - l^ijV ■ u, (32) 

where commas denote partial differentiation. Note that the form of the stress tensor 
above applies only to a monatomic gas. In more general cases there may be additional 
contributions from the bulk viscosity corresponding to terms proportional to (^jV • u. 
To compare with the incompressible case, we evaluate 

- V • r = v \v 2 u + \W ■ u + 2S • V ln(pu)] , (33) 
P J 

and note that, in addition to the V 2 u term there is also a term VV • u, which vanishes 

in the incompressible case, and a term Sy V 3 - ln(pu), which vanishes when the dynamical 

viscosity, p = pu, is constant. 



Equation (31) has to be solved together with the continuity equation 



^ = -pV. U , (34) 

and an energy or entropy equation, 

pT^ = 2upS 2 . (35) 

The heating term is generally given by u^jTij. Splitting Uij = + a^- into symmetric 
and antisymmetric parts, it is clear that only contributes after multiplying with 
another symmetric matrix, i.e. with ry. Furthermore, since 7y is also trace- free, the 
result does not change when adding or subtracting from a term proportional to 5^, 
in particular ■ u. Therefore, we have u^jTij = 2pz/S 2 , which is manifestly positive 

definite. 

For a perfect gas the specific entropy s is related to pressure and density via 

s = c v hip — c p liap + s , (36) 

where s is an additive constant. (The specific entropy s is not to be confused with 
or Sy.) It is important to realize that even in the inviscid limit, v — > 0, the term 
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2upS 2 cannot be neglected in equation (35). For example across a shock there is always 
a well-defined increase in specific entropy that is independent of the value of v. 

In compressible fluid dynamics it is often advantageous to consider the evolution 
equations in their conservative form. This means that the rate of change of the density 
of a conserved quantity, X, is given by the negative divergence of its corresponding flux, 
i.e. 

d 

— (density of X) = —V • (flux density of X) + sources — sinks, (37) 

L/ v 

where the presence of sources and sinks indicates additional processes whose detailed 
evolution is not captured by the total energy equation within the same framework. An 
example is radiation, which provides sources and sinks to the energy equation as heating 
and cooling terms. Alternatively, if the evolution of the radiation energy is included in 
the total energy equation, any explicit heating and cooling terms disappear, and only 
boundary (flux divergence) terms remain |114] . If there is no radiation, gravity, external 
forcing, etc, there are no additional terms, so the conservative form of the equations is 

i = -<| ( ^ < 38) 

d d 

— {pUi) = -— {pu.Uj + dijp - Tij) , (39) 
d d 

— (pe + \pu 2 ) = -— (pujh + \pUjU 2 - u^) , (40) 

where h = e + p/p is the specific enthalpy per unit mass. For a perfect gas, h and e are 
proportional to temperature with h = c p T and e = c v T, where c p and c v are the specific 
heats at constant pressure and constant volume, respectively. 

The equations above show explicitly that the volume integrals of the terms under 
the time derivative are conserved, i.e. constant in the absence of fluxes in or out of the 
domain. In one dimension, the terms in parentheses under the spatial derivatives are 
constant and, in particular, uniform across a shock. This allows shock jump conditions to 
be derived. Note that, since viscosity acts only locally, these conditions are independent 
of the width of the shock. This is an important property that allows simulating 
highly supersonic turbulence using a modified viscosity (Neumann-Richtmyer artificial 
viscosity) for smearing out the shock |115j . In the presence of source or sink terms in 
equations (38)-(40) this would no longer be possible. 



4-3. Anelastic approximation 

The advantage of making the assumption of incompressibility is not only that one has 
one equation less to solve (the dp/dt equation), but mainly that one eliminates sound 
waves, whose associated wave speed is often much faster than the speed associated with 
other processes. This means that one can then focus more efficiently on the slower 
dynamics of the system. 
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Incompressibility is normally associated with constant density. In view of our earlier 
discussion regarding the strong density stratification in stars, incompressibility would 
not be a useful assumption, even though the sound speed can be much larger than other 
speeds such as that associated with the convection itself. It is then better to relax the 
condition = and use instead V-pw = 0. This is called the anelastic approximation 
|116[ 1117] . It is important to realize that with this assumption one replaces the original 



continuity equation (38). Consequently this equation can then no longer be used to 
argue that dp/dt = 0. Indeed, p is in general not constant in time and can evolve, while 
V • pu = is maintained at all times. This technique is sometimes used in simulations 
of solar convection [TTBl EES H2H ESI H2D H22| - 

Just like in the incompressible case, also here one has to solve a Poisson-like equation 
that emerges when taking the divergence of the evolution equation for the momentum 



density. Taking the divergence of equation (39) one obtains 

V 2 p = V ■ R, (41) 
where R = —pu ■ Vm + F + V • r is the sum of the advection term plus all the other 



terms on the right hand side of equation (39), except for the pressure gradient term 



The F term in equation (41 ) refers to additional terms such as gravity and Lorentz force 



terms in equation (31). 



The anelastic approximation is sometimes associated with linearizing the equation 
of state [55J. However, this is not necessary and one can just continue working with the 
original, fully nonlinear equation of state |119[ 1123j . The only difference is that in the 
fully compressible case one would obtain the pressure from density and specific entropy, 
while in the anelastic case one obtains the density from pressure and specific entropy, if 
the latter is indeed the main thermodynamic variable. 

4-4- Large eddy and hyperviscous simulations 

The maximum achievable Reynolds number scales as the number of mesh points in 



one direction, raised to the power 4/3; see equation (11). With the largest attainable 
resolution being at present 4096 3 [18J, it is impossible to reach Reynolds numbers of 
10 6 and beyond. In many engineering applications of turbulence one needs to calculate 
flows at very large Reynolds numbers and one therefore uses large eddy simulations. 
This involves some representation of the unresolved Reynolds stress in terms of other 
flow variables. This approach can be rather uncertain. Unlike engineering applications, 
where such models can be tested against measurements, this is usually not possible in 
astrophysics, due to a large number of additional complications (strong stratification, 
magnetic fields, rotation, etc.) that are hard to realize in the laboratory. The best one 
can therefore hope for is a rigorous comparison of large eddy simulations with DNS. 
Examples of this are discussed in §[6j 

One of the simplest subgrid scale models is the Smagorinsky model |124j . This 
approach is strictly dissipative, i.e. the Reynolds stress of the unresolved velocity 
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fluctuations, denoted here by primes, is modeled by a viscous stress of the form 



uty = -2(C s Ax) 2 p|S|S ii , (42) 

where Cg is the Smagorinsky constant (between 0.1 and 0.2) \12h\ I126j 1127] . and the 
rate-of-strain tensor S was defined in equation ( 32 ) , and is here applied to the resolved 
motions u, i.e. excluding the subgrid scale motions. Another approach, which cannot be 
classified as large eddy simulation, consists in using hyperviscosity. In spectral space, 
the viscosity operator —vk 2 is simply replaced by —u n k 2n , where n > 1 is the order 
of hyperviscosity. Unlike the Smagorinsky model, the results from this approach are 
known not to converge to the original Navier-Stokes equations, but the hope is that in 
the inertial range the flow is unaffected by the unphysical form of the diffusion operator. 
This is indeed the case, as was demonstrated in Ref. [104] . 



4-5. Turbulence simulations using Godunov/PPM-type schemes 

The Godunov scheme is a conservative numerical scheme for solving partial differential 
equations. In this method, the conservative variables are considered as piecewise 
constant over the mesh cells at each time step and the time evolution is determined 
by the exact solution of the Riemann shock tube problem at the intercell boundaries. 
This scheme consists of first defining a piecewise constant approximation of the solution 
at the next time step. The resulting scheme is usually first-order accurate in space. 
This approximation corresponds to a finite volume method representation whereby the 
discrete values represent averages of the state variables over the cells. Exact relations for 
the averaged cell values can be obtained from the integral conservation laws. Next, the 
solution for the local Riemann problem is obtained at the cell interfaces. This is the only 
physical step of the whole procedure. The discontinuities at the interfaces are resolved 
as a superposition of waves satisfying locally the conservation equations. The original 
Godunov method is based upon the exact solution of Riemann problems. However, 
approximate solutions can be applied as an alternative. Finally, the state variables are 
averaged after one time step. The state variables obtained after the second step are 
averaged over each cell defining a new piecewise constant approximation resulting from 
the wave propagation during the time step. 

Nowadays one uses often higher-order Godunov schemes for astrophysical 
applications. One such method is the piecewise parabolic method that is also referred 
to as PPM. Examples of such codes include Athena |128j . Pluto |129j . Nirvana |130j . 
Ramses |131j . Flash |132j . and Enzo |133j . Such codes have been used for many 
astrophysical applications including supersonic, isotropic homogeneous turbulence |134] . 



4-6. Analyzing and modeling turbulence with wavelets 

Wavelets are sometimes used both to analyze and to model turbulence. In particular 
the wavelet technique has been used for extracting coherent vortices out of turbulent 
flows. The aim is to retain only the essential degrees of freedom responsible for the 
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transport. It is intriguing that with this technique one can actually retain nearly all 
velocity structure and dissipation information in turbulent flows by using a relatively 
small selection of wavelets with non-zero amplitudes |135j ; see also Ref. |136j . This 
method is related to the so-called proper orthogonal decomposition of turbulent flows 
[137J. This decomposition is statistically based and permits the extraction of spatio- 
temporal structures that are judged essential according to predetermined criteria. It is 
not only useful in the analysis and synthesis of data from simulations and experiments, 
but it also allows the construction of low-order models from the Navier-Stokes equations. 
Finally, let us note that the wavelet representation has been applied with success to 
simulations of resistive drift- wave turbulence in magnetized plasma Hasegawa-Wakatani 
system |138j . 



5. Extra ingredients to turbulence in astrophysical flows 

5.1. Passive scalars: mixing and dust dynamics 

One of the simplest additional ingredients in fluid dynamics in general, and in turbulence 
physics in particular, are passive scalars. The passive scalar concentration per unit mass, 
9, is governed by the equation 

d d ( d9 \ 

(pO) = - 7 r-[pu ] 9-p Ke —), (43) 



dt dxj \ dx d f 

where kq is a diffusion coefficient for the passive scalar concentration. This equation 
describes the transport of chemicals in a gas. Additional source and sink terms could 
be included to model production and destruction of chemicals. The non- conservative 
form of this equation can be written as 

Y)9 1 

— = —V ■ (p Ke V9) , (44) 



where we have made use of the continuity equation (34). For Kg = 0, this equation gives 
D9/T)t = 0, which shows that the concentration per unit mass is unchanged at each 
point comoving with the flow. 

Another class of scalars are inertial particles that are advected by their own velocity 
u p rather than the velocity of the gas u. The evolution equation of u p is similar to that 
of u, except that it lacks the pressure gradient term and the Lorentz force. However, 
such particles are strictly speaking active particles, because of the mutual coupling 
between the two velocity fields. Only in the limit of sufficiently light particles can the 
back-reaction on u be neglected. 

In astrophysics one often finds the condensation of heavier elements into solid dust. 
Their evolution is described as a passive scalar or as passively advected particles. The 
inclusion of inertia can sometimes become important, because inertial particles have a 
tendency to accumulate in anti-cyclonic vortices [1391 11401 11411 11421 1143] . 
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5.2. Active scalars: stratification and convection 

In this context, the term "active" refers to the property that the scalar quantity can 
affect the momentum equation, for example by exerting a pressure gradient force. An 
example is the advection-diffusion equation for the energy density of low-energetic 
cosmic rays |144l [145J. Another example concerns temperature or specific entropy, 
which affect the momentum equation by locally changing the relation between pressure 
and density. In the presence of a gravity force, F = pg, this can lead to an Archimedian 
buoyancy force. Furthermore, with g ^ a new wave mode can exist known as gravity 
waves (not to be confused with gravitational waves of the space-time metric in general 



relativity; see comment at the end of § 3.7). The restoring force comes from the linearized 
buoyancy term, (Sp/p )g ~ (Sp/p + 5s/c p )g. Since the restoring force is related to 
gravity, these wave modes are often referred to as g-modes, in contrast to p-modes 
or sound waves, whose restoring force is related to the pressure gradient. If pressure 
fluctuations may be neglected the essential terms are 

-^ = ... + 5sg/c p , (45) 
— = -u — (46) 

dt - Uz dz' 1 j 

where s denotes the specific entropy of the background stratification. The oscillation 
frequency Abv (for Brunt- Vaisala frequency) is given by 

A^bv = -9 ■ Vs/c p . (47) 

While this pair of equations represents the basic feedback loop correctly, it ignores 
the fact that buoyancy is only possible when there is lateral non-uniformity of density. 
Indeed, solving the proper dispersion relation reveals that on large scales the frequency 
increases linearly with wavenumber; see, e.g., Ref. |146] for a review. In Figure[8]we show 
the dispersion relation as a function of the horizontal wavenumber, k r = (k 2 + k 2 ) 1 ^ 2 , 
for k z = and different values of the ratio of specific heats ranging from 7 = 1.1 to 
1.9. The p-modes correspond to the upper branch while the g-modes to the lower one. 
Also shown are the g-modes obtained using the anelastic approximation discussed in 



4.3| Note that this approximation yields correct results for 7 close to one and for large 
horizontal wavenumbers, i.e. on scales that are small compared with the pressure scale 
height HH1. 

Given that gravity points downward, is positive (i.e. the frequency is real) 
when the specific entropy increases in the upward direction. If it decreases with height, 
the system is unstable to the onset of convection with an approximate growth rate 
given by Im|A^Bv|- Here we have omitted viscous and diffusive effects that could slow 
down the growth and even stabilize the system. This is quantified by the value of the 



Rayleigh number that will be defined and discussed in more detail in § 6/7 However, 
in astrophysics viscosity and diffusivity are comparatively small and one uses just the 
condition g ■ Vs > for instability. This is known as the Schwarzschild criterion and 
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Figure 8. Dispersion relation uj = u)(k r ), for k z = with k r = (k* + A:^) 1 / 2 
being the horizontal wavenumber, for different values of the ratio of specific heats, 
7, showing p-modes branch (upper branch) and g-models (lower branch) compared 
with the case where the anelastic approximation has been made (dashed line) . Length 
is given in units of Hq = ^fH p /(l — 7/2) and time in units of To = Hq/c s , and 
ujq = y/ry — 1/(1 — 7/2) is the non-dimensional Brunt -Vaisala, frequency. In the 
plot, the break point at k r /uo = 1 corresponds to a critical horizontal wavelength 



of 4 r it = 27T7-ff p /V7 - 1- 
Ref. [137] . 



For 7 = 5/3 this means £ C rit ~ 12.8 i? p . Adapted from 



corresponds to saying that the Rayleigh number is positive (convection is discussed in 



more detail in §6.7). 



In the presence of strong vertical density stratification, the convection flow tends to 
develop an interconnected network of downdraft lanes, with isolated tube-like stronger 
downdrafts at network vertices. With depth, the downdrafts merge and the network size 
increases [148] ; cf. also Figs. |3]and|4j At large Reynolds numbers the flow is of course 
turbulent, but with the intensity of turbulence strongly influenced by stratification 
effects: Because ascending flows are strongly divergent, turbulence is suppressed there, 
while in downflows, which are converging, turbulent intensity is enhanced. 

In the Sun, the Prandtl number, Pr = u/x, is far below unity (around 10 -5 ). 
This means that velocity or vorticity structures can be much thinner than temperature 
structures. As a consequence, thin vortex tubes can develop within downdrafts. The 
dynamical pressure associated with vortex tubes allows locally a lower gas pressure and 
hence a lower density, making vortex tubes buoyant. As a result, the downdraft speed is 
slowed down ( "vortex braking" ) [1491 1150] . This is a particular property of low Prandtl 
number dynamics which, at the same time, requires compressibility. 
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Compressibility leads to yet another interesting effect in convection. The pressure 
gradient associated with driving the horizontal expansion of upwelling motions works 
in all directions, and in particular also in the downward direction. This tends to brake 
the upwellings. This phenomenon is known as buoyancy braking [151 J. 

Another important effect caused by compressibility is the production of vorticity by 
the baroclinic term, i.e. the curl of p _1 Vp. The curl of this term is finite if the surfaces of 
constant p and p are inclined relative to each other. Another way of writing this term is 
by using the thermodynamic relation for the differential of enthalpy, dH = TdS + Vdp. 
With this we can write the pressure gradient term in terms of specific enthalpy, specific 
entropy s, and specific volume p -1 as 

- p- l Vp = -Vh + TVs. (48) 

This formula shows that the baroclinic term is just given by 

V x (-p^Vp) = VT x Vs. (49) 

This relation will becomes useful later in connection with the Taylor-Proudman theorem 
and ideas to understanding departures from it. The baroclinic term vanishes under 
isothermal (T = const), isentropic (s = const), or barotropic [p = p(p)] conditions. In 



all these cases, equation (48) can be written purely as a gradient term, — Vh, where h is 
then called the pseudo-enthalpy and it is proportional to h which, in turn, is proportional 
to the temperature. In the irrotational case, uj = 0, the only nonlinearity comes from 
the \u 2 term in the reduced pressure. 

5.3. Rotation and shear 

It is often convenient to solve the governing equations in a rotating frame of reference. 
In that case, Coriolis and centrifugal forces as well as possibly the Poincare force have 
to be included on the right hand side of the Navier-Stokes equation, so the equation 
takes the form 

Du 

= ■■■ — 2f2 x u — f2 x (f2 x r) — f2 x r, (50) 

where r is the position vector with respect to a point on the rotation axis and Qq = const 
is the angular velocity vector of the reference frame. The Poincare force, Ct x r can 
drive flows and even turbulence in precessing bodies with boundaries. This has been 
discussed in attempts to explain the flows that drive the geodynamo |152[ I153[ I154[ 1155] . 

An important effect of the Coriolis force is to suppress variations of the azimuthal 
velocity in the axial direction. This can be seen by taking the curl of the Coriolis term, 

V x (-2Q x u) = 2fi ~ £V ' ' ( 51 ) 

where u± — u — (u ■ z)z is the velocity in the direction perpendicular to the rotation 
axis and z is the unit vector along the direction of f2 . Taking the curl of the evolution 
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equation for du/dt in cylindrical polar coordinates, (R,4>,z), and projecting on the <fi 
component, yields 

* ' w = 2fio §? " * ' (VT x Vs) + (52) 

where FJ 1SC is the component of the viscous force. This shows that, when rotation is 
important, i.e. fl is large, du^/dz must be small. 

Of course, the physics is independent of the coordinate system one is working in. 
If one is working in the inertial (non-rotating) frame, there is no Coriolis force, but 
one can then write = RQ, where Q = Q(R, <p, z) is now the local angular velocity, 
which is not to be confused with f2o- If the velocity has only an azimuthal component, 
u = (pRQ, one can write the curl of the u ■ Vm term as 

dVL 2 

V x (-u-Vu) = R—. (53) 

oz 



We will return to the astrophysical consequences of this in §7^3 when we discuss the 
angular velocity of the Sun. 

5.4- Active vectors: magnetic fields and dynamos 

An important vector field to be included in the fluid equations is the magnetic field, 
B. It is an active vector because the Lorentz force, J x B, backreacts through the 
momentum equation, so 

p^ = ... + JxB, (54) 

where J is the current density. One makes here the assumption that there is no net 
charge in the fluid, i.e. the density of positive and negative charge carriers balances 
everywhere, and the currents are produced by the sum of the fluxes of counterflowing 
positive and negative charge carriers. The B field is solenoidal and its evolution is 
governed by the Faraday equation, 
dB 

— = -V x E, with V • B = 0, (55) 
where E is given by Ohm's law, 

- E = ux B - J/a, (56) 

and a is the electric conductivity. Its inverse is related to the magnetic diffusivity, 
7] = (/io°") _1 , and it has the same dimension as the kinematic viscosity v. Their ratio is 
the magnetic Prandtl number Pr^ = v jr\. 

Ampere's equation is used to express the current density in terms of the magnetic 
field via 

J = V x B/tJ, Q , (57) 



where /iq is the vacuum permeability. Equation (57) is an approximation to the 



full Faraday equation which includes also the displacement current. Neglecting it 
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A-BdV = -2a~ 1 / J-BdV- fF H -dS, (60) 



corresponds to filtering out electromagnetic waves, which is justified at finite electric 
conductivity and velocities small compared with the speed of light. 

Inserting equation (56) into equation (55) we obtain the induction equation in the 

form 

f) T-t 

— = Vx(uxB- J I a) . (58) 

In its 'uncurled' form this equation reads 
dA 

— = ux B - J jo - V0, (59) 

where <fi is the electrostatic potential. By evaluating the time derivative of A ■ B and 
integrating over space we obtain the evolution equation for magnetic helicity, 

d 

eft 

where Fjj = E x A + <pB is the flux of magnetic helicity. 

Magnetic fields constitute an additional form of energy, Em = J B 2 /(2fio) dV, 
whose evolution is given by 

dll 27T dV = ~ J U '( J X B ) dV j J 2 dV - <f F M -dS, (61) 

where F M = E x B / fi is the Poynting flux. Equation (40) for the evolution of the 
total energy density can be generalized correspondingly by adding B 2 /2fi underneath 
the time derivative and Fm underneath the divergence term. 

In this connection it might be useful to emphasize that in numerical simulations 
one hardly uses the full energy equation in that form if the magnetic energy becomes 
comparable to or in excess of the thermal energy. Normally one would calculate the 
thermal pressure from the internal energy, but in the magnetically dominated case this 
becomes a small residual between total, kinetic, and magnetic energies, and so this 
calculation becomes exceedingly inaccurate. 

Another comment regarding simulations is here in order. A commonly encountered 
difficulty is to preserve solenoidality of B. One method is to use a staggered mesh 
and to evaluate the right-hand-side of equation (55) such that the numerical evaluation 
of the curl produces zero divergence to machine accuracy. Another method is to use 
A as dependent variable, which also preserves V • B = 0, and it also allows for a 
straightforward calculation of the magnetic helicity. Yet another method is to write 

B = Vax V/3, (62) 

where a and /3 are the Euler potentials |156j . However, this method only works in the 
strictly ideal case, in which case the evolution equations are just 

Da/Dt = D/3/Dt = 0. (63) 

This approach is now quite popular in smoothed particle hydrodynamics calculations, 
because then the values of a and are just kept fixed at each Lagrangian particle 
|157l 1158] . Unfortunately, this method cannot even approximately capture non-ideal 
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effects. As a consequence, dynamo action (see below) is not possible in this approach 
and energy spectra of MHD turbulence with imposed field become too shallow |159j . 
Finally, there is the possibility of divergence cleaning, which requires the solution of a 
Poisson-type equation for the correction term to the numerically obtained B field. This 
approach is analogous to calculating the pressure under the constraint that V • u = 
or V ■ pu = 0; see §4^3 The disadvantage here is that this approach may introduce an 
unphysical nonlocality as a consequence of invoking a Poisson-type equation. 

The Lorentz force gives rise to various restoring forces that lead to additional wave 
forms including Alfven waves as well as fast and slow magnetosonic waves. The slow 
magnetosonic waves are particularly important in the presence of shear and rotation, 
because those waves can be destabilized to give rise to the magneto-rotational instability. 



This will be discussed in more detail in 8 6.11 



One of the other new features allowed by the addition of magnetic fields is the 
possibility of self-excited dynamo action, i.e. the spontaneous conversion of kinetic 
energy into magnetic energy by work done against the Lorentz force. This is an 
important process in astrophysics. Magnetic fields observed in planets and stars with 
outer convection zones are clear examples where dynamo action is required to sustain 
magnetic fields against ohmic decay and to explain field reversals on time scales short 
compared with the resistive time. Galaxies and clusters of galaxies also harbor magnetic 
fields. Many spiral galaxies show magnetic fields with a large-scale design that is 
approximately axisymmetric. One prominent exception is a galaxy with the name M81, 
where the field is non-axisymmetric with a strong m = 1 component, i.e. the field is 
proportional to e 1Tn ^, where <fi is the azimuthal angle. Observations give direct indications 
about the turbulent nature of galactic discs, so the magnetic field must be maintained 
against turbulent decay in the vertical direction along the axis. The relevant time scale 
is only about 10 7 yr. In the present review we discuss dynamos only insofar as they are 
directly connected with understanding or clarifying astrophysical turbulence. 

Details regarding dynamo theory as well as magnetic fields in solar-like stars and 
galaxies have recently been reviewed in Refs. [TTI I5B* | 1160] . One of the important recent 
developments concerns the realization that the evolution of the large-scale magnetic field 
can be constrained decisively by magnetic helicity evolution; see equation (60). This 
has to do with the fact that large-scale magnetic fields tend to be helical. This point 



will be taken up briefly in § 7.2, but for a more thorough discussion we refer to Ref. [11] 
for a recent review. 

In the incompressible case with constant density p = po, it is convenient to write the 
MHD equations using Elsasser variables z± = u±B / \/poPo, because then the evolution 
equations take a form similar to the usual Navier-Stokes equations, i.e. 

^ + Zt ■ Vz± = -VII + zA7 2 z±, V • z± = 0. (64) 
at 

Here, v = rj has be assumed for simplicity, and n = (p + B 2 /2p )/p is a pressure that 
ensures that V • z± = 0. 
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5.5. Radiation: optically thick and thin 

Radiation transport describes the coupling to the photon field. As far as the dynamics 
is concerned, the radiative flux gives rise to a radiation force that can for example 
cause levitation of the gas by radiation. The radiative energy flux divergence enters 
the energy equation and describes local heating and cooling. Thus, the momentum and 
specific entropy equations are amended as follows, 

T)u on „ . 

"m = ■ + IT dl (65) 

pT^ = ... - V ■ F„ d . (66) 

Here, k is the opacity, i.e. the photon cross-section per unit mass. The cross-section per 
unit volume is pn, which is also the inverse mean free path of photons, i = (pn)^ 1 . If 
the mean free path is small compared with other relevant length scales, a diffusion 
approximation may be used for -Frad, which means that it is proportional to the 
negative gradient of the radiation energy density, F r!3id = — |£cV(aT 4 ), and so it points 
in the direction of the negative temperature gradient. The transition layer between 
optically thin and thick is an important region in astrophysics, because it marks the 
effective surface of an otherwise extended body. In this transition region the diffusion 
approximation is no longer valid and proper equations for the radiation intensity have 
to be solved to obtain .F rad ; see Refs. |119^ I161j . 



6. Simulations of turbulence 



Astrophysical turbulence is frequently caused by instabilities. However, many 
instabilities imply the presence of anisotropies. For example in convection the vertical 
direction is a preferred one, while in the case of the magneto-rotational instability 
the velocity gradient matrix associated with the shear governs the anisotropy. In the 
presence of magnetic fields, the otherwise isotropic turbulence becomes at least locally 
anisotropic, because at every patch in the turbulent flow the direction of the local mean 
field imprints anisotropy on all smaller scales within this patch. On the other hand, 
much of turbulence theory is concerned with isotropic turbulence. Computationally, 
isotropic turbulence can be modeled by adopting an imposed forcing function. Common 
applications of isotropically forced turbulence include simulations of turbulent star 
formation, as well as turbulent mixing and dynamo processes. We begin by discussing 
some general aspects of isotropic turbulence simulations. 



6.1. General aspects 

The concept of isotropic turbulence is a convenient and useful theoretical idealization. 
Computationally, isotropy does not lead to any significant simplification, except that 
periodic boundary conditions are possible and in many ways advantageous. Isotropic 
turbulence needs to be forced by an isotropic body force, unless an isotropic instability 
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can be identified that would drive turbulence. The thermal instability would be an 
example of an instability without preferred direction, but simulations have not shown 
that it can lead to sustained turbulence |162[ I163j . The Jeans instability is another 
example, which is particularly relevant to the problem of star formation through strong 
compressions by the turbulence in the interstellar medium. This problem is frequently 
tackled using smoothed particle hydrodynamics |164[ 1165] , while mesh-based techniques 
have explored mostly the case of forced supersonic turbulence |166} 1167] and have only 
recently incorporated the effects of self-gravity, augmented with so-called sink particles 
to account for the production of high density cores that cannot be resolved with a fixed 

mesh pansm 11701 nrg . 

In order to study more basic properties of turbulence one often resorts to a random 
forcing function to simulate the effects of an instability with a well-defined forcing 
strength and a well-defined length scale of the driving. Plane waves with randomly 
changing orientation is an obvious possibility for driving turbulence. To make the forcing 
divergence-free, one uses only transversal waves. 

The idea of simulating turbulence on the computer developed during the 1970ies. 
Almost all simulations in those days utilized pseudo-spectral methods, i.e. spatial 
derivatives are calculated in Fourier space by multiplication with ik, but all nonlinear 
terms are calculated in real space. The main advantage of such methods is the small 
discretization error. Furthermore, this technique also allows an efficient solution of the 
Poisson-like equation for the pressure if one makes the incompressible or the anelastic 
approximation, i.e. V • u = or V • pu = 0, respectively. 

Spectral methods have the disadvantage that one cannot easily deal with arbitrary 
boundary conditions. Also, the Fourier transformation is a nonlocal operation which 
is not optimal when using many processors. These are reasons why sometimes finite 
difference methods are used instead. Finite difference methods are normally not as 
accurate as spectral methods unless one uses a higher order scheme (e.g. fourth and 
sixth order schemes are common choices). On the other hand, many astrophysical flows 
develop shocks for which there are a number of other dedicated methods (Riemann 
solvers, approximate Riemann solvers, monotonicity schemes, Godunov schemes, and 
Neumann-Richtmyer artificial viscosities |115] 1172] ). These methods are frequently 
generalized to mesh refinement methods that allow increased accuracy in specific 
locations of the flow. Finally, there are also Lagrangian methods of which Smooth 
Particle Hydrodynamics is an example [1651 11731 11741 1175] . A promising new Lagrangian 
method has been presented in Ref. |176j . 

6.2. Hydrodynamic turbulence 

When simulations became able to resolve turbulence with around 128 3 meshpoints, it 
became evident that much of the flow is governed by a tangle of vortices; see, e.g., 
Refs. |177] I178[ 1179] . The left-hand panel of Figure [9] shows examples of such vortices. 
Their thickness is related to the viscous scale while their length was often expected 
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Figure 9. Examples of vortex tubes in homogeneous turbulence. Courtesy of Zhen-Su 
She (left figure) [177] and Paul Woodward (right figure) fL79] . 



to be comparable with the integral scale. However, in subsequent years simulations at 
increasingly higher Reynolds numbers seem to reveal that the vortex turbulence become 
a less prominent feature of otherwise nebulous looking structures of variable density (see 
the right-hand panel of Figure [9]) . 

Incompressible forced turbulence simulations have been carried out at resolutions up 
to 4096 3 meshpoints [IB] . Surprising results from this work include a strong bottleneck 
effect [20] near the dissipative subrange, and possibly a strong inertial range correction 
of about A; -0 ' 1 to the usual A; -5 / 3 inertial range spectrum, so that the spectrum is A;~ L77 . 
Note that the She-Leveque correction (17) is only A; -0 03 , so that the spectrum is A;" 1,70 . 
Similarly strong inertial range corrections have also been seen in simulations using a 
Smagorinsky subgrid scale model [19] (512 3 meshpoints, dashed line in Figure 10). 
Here we also show the results of simulations with hyperviscosity, i.e. the z/V 2 diffusion 
operator has been replaced by a ^V 6 operator |180] (512 3 meshpoints, dash-dotted 
line). Hyperviscosity greatly exaggerates the bottleneck effect, but it does not seem to 
affect the inertial range significantly; see Figure [TQ| 

The nature of the k~ 0A correction factor is currently not understood. It might be 
an artifact resulting from applying a forcing function at a scale close to the scale of the 
box [S\. Alternatively, the presence of a bottleneck might also lead to the emergence of 
a dip just before the bottleneck. In either case this would not be a true k~ 0A correction 
in the entire inertial range. 

In virtually all astrophysical settings the relevant Reynolds numbers are very large 
and the bottleneck is hardly important, because it is located at very small length scales. 
However, this is not the case in simulations which show the bottleneck as a pronounced 
feature. There are several important issues here. Firstly, simulations at resolutions of 
256 3 meshpoints give hardly any indication of a bottleneck effect, and only at resolutions 
of 1024 3 meshpoints and above does it really develop its full strength. For this reason 
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Figure 10. Comparison of energy spectra of the 4096 3 meshpoints run [IB] (solid line) 
and 512 3 meshpoints runs with hyperviscosity (dash-dotted line) and Smagorinsky 
viscosity (dashed line). (In the hyperviscous simulation we use 1/3 = 5 x 10~ 13 .) 
The Taylor microscale Reynolds number of the Kaneda simulation is 1201, while the 
hyperviscous simulation of Ref. [180 has an approximate Taylor microscale Reynolds 
number of 340 < Re^ < 730. For the Smagorinsky simulation the value of Re,*, is 
slightly smaller. Courtesy of Nils E. Haugen 19J. 



the bottleneck effect has been studied more seriously only in recent years. Secondly, 
the bottleneck effect can affect certain aspects of a simulation in a way that is not 
yet asymptotically meaningful. An example is the small-scale dynamo effect that is 
discussed below. 

6.3. Supersonic turbulence 

In the interstellar medium the gas can condense into more concentrated regions called 
molecular clouds. These clouds are so cold that molecules can form, which explains their 
name. Because of low temperature in the range of 10-100 K, the flows in these clouds can 
become highly supersonic. This in turn leads to even stronger mass concentrations that 
can become gravitationally unstable and form stars. This is why supersonic turbulence 
is commonly studied in connection with star formation [1671 1131 j . 

With increasing Mach number, density fluctuations begin to become important. 
In fact, in supersonic turbulence with an isothermal equation of state it has been 
demonstrated that the standard deviation of the (linear) density, o"ii nea r, grows linearly 
with the Mach number [Ml HH21 UBS] 

^linear = 7 Ma > (67) 

where the Mach number is defined as Ma = u rms /c s . The density obeys a log-normal 
distribution, i.e. the probability density function, p(lnp), with Jp(lnp)dlnp = 1, is 
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Figure 11. Left: time- averaged velocity spectra compensated by k 2 from uniform 
grid PPM simulations at resolutions 256 3 , 512 3 , 1024 3 , and from an adaptive mesh 
refinement simulation with effective resolution of 2048 3 grid points. The spectra 
demonstrate convergence for the inertial range of scales. Right: time-averaged 
spectrum of p x l z u compensated by fc 5 / 3 . The straight lines represent the least-squares 
fits to the data within a suitable intermediate range of wavenumbers. Courtesy of 
Alexei Kritsuk [15] , 



given by 



p(lnp) 



1 



exp 



|(lnp-lnp) /a 5 



(68) 



V2710- 2 

where a is the standard deviation of the logarithmic density, which is related to the 
Mach number like [166J 



a 



again with 7 



= In (l + 7 2 Ma 2 J , 
1/2 to good accuracy. 



(69) 



As indicated in §2.1 the spectra of u, p l l 2 u, and p l / 3 u begin to differ from each 
other at larger Mach number. Observations of the line-of-sight velocity dispersion of 
molecular clouds in the Perseus cluster also show that in the highly supersonic case 
the velocity spectrum is not far from k^ 1 - 8 [184], and thus deviates clearly from the 
characteristic spectrum of shock turbulence |185] . However, the density weighted spectra 
tend to become shallower. In particular the spectra of p 1 ^ 3 u are very close to k~ 5/3 P]; 



see Figure VL This appears to be connected with the fact that the kinetic energy flux, 
i.e. the quantity that is constant throughout the inertial range at scale I, is given by 
puf/l. This idea goes back to an early paper by Lighthill |12j . 



6.4- Hydro-magnetic turbulence 

The gas in many astrophysical settings is partially or fully ionized and hence electrically 
conducting. This means that the effects of magnetic fields cannot be neglected. 
The full extent of associated behaviors is not yet well understood, nor is there 
unambiguous evidence for universal and asymptotic scaling behavior in the limit of 
large fluid and magnetic Reynolds numbers [186J. However, using decay simulations at 
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moderate magnetic Reynolds numbers, three types of behavior have been identified |187] , 
depending essentially on the ratio of the initial magnetic to kinetic energy densities. The 
purpose of this section is to discuss the expected properties in these three regimes. 

It is convenient to introduce here the Alfven speed va = B rms / ^//ioPo associated 
with the random magnetic field -B rms . The case of sub-equipartition random fields with 
< Mrms was studied by Iroshnikov |188] and Kraichnan [189 J who argued that the 
turbulence can still be treated as isotropic and that the flux of energy e down the 
turbulent cascade will be modified by Alfvenic interactions and replaced by the geometric 
mean of energy flux and Alfven speed, i.e. 

e {ev A ) 1/2 . (70) 

The dimensional argument used in equation ^ for the energy spectrum of Kolmogorov 
turbulence gets correspondingly modified and is then of the form 

E(k) = C IK (ev A y/ 2 k-V 2 . (71) 

In the case of strong magnetic fields, va 3> Mrms, the turbulence becomes highly 
anisotropic, so the spectrum E(k±,k\\) depends on the wavenumbers parallel (k\\) and 
perpendicular (k±) to the local direction of the field. In this limit the turbulence can 
be treated as wave turbulence using weak turbulence theory |190j . which leads to 

E{k ± ,k\\)~kZ 2 . (72) 

In the intermediate case, the kinetic energy of the turbulence is comparable to that of 
the magnetic field. This regime is referred to as strong turbulence - not because the 
field is strong, but because the u ■ Vm nonlinearity cannot be neglected. The flow is still 
anisotropic, and energy is cascaded in k± at a rate e. The resulting energy spectrum is 
[mi 11521 [1951 1153] . i.e. 

£(fc±,fc||) = C GS e 2 /V /3 . (73) 

In the following we present a more detailed phenomenology that highlights the essential 
physics behind the various regimes. 

We consider as governing equations the MHD equations written for the Elsasser 
variables z±, see equation (64), and denote by the modulus of z± at wavenumber 
k±. In all cases the energy spectrum is given by 

E{k ± ,k\\)~zljk ± , (74) 

and the spectral energy flux is then given by an expression of the form 

e = Z k ± / T casc, (75) 

where r casc is the cascade time. The main difference between the various regimes lies in 
the form of the r casc ; see also Refs. [HI I195j . 

For strong magnetic fields, interactions are being accomplished by wave packets 
traveling in opposite directions. The duration of the interactions is given by 

T A ={vAh)-\ (76) 
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where fen 1 is the longitudinal extent of such a packet. The fractional change in a wave 
packet is given by the ratio 

X = W t nl (77) 
of Alfven time to the nonlinear interaction time 

tnl = {z k± k x )~ l - (78) 
However, because the sign of each interaction is random, the net effect grows only like 
the square root of the number of interactions. Therefore, the effective fractional change 
associated with each interaction is only \ 2 - This means that the effective cascade time 
is T casc = ta/x 2 - By contrast, in the strong turbulence regime the Alfven and nonlinear 
times are equal, i.e. \ — 1? an d the cascade time is therefore just r casc = r NL - Since 
the zu ± in equation (74) enters also in the expression for r casc , the resulting spectra are 
qualitatively different. For weak turbulence we have 

e = ^ = zl i ^^, so Eik^-^-ieu^WkZ*, (79) 

Tcasc ^A^j K± 

while for strong turbulence we have 



2 2 

k 



Z k ± 2 - - - 



-7. (z k± k±), so £(fe ± ,fe|,)~-^~e 2 / 3 fe7_ 5/3 . (80) 



In the latter case, because of r NL = ta, we have k±/ku = va/ Zk ± = (vx/ejk 1 / 3 , so 
the degree of anisotropy increases toward smaller scales until we have fen — > oo. For 
weak turbulence we have fen — > 0, so the turbulence is fully anisotropic at all scales. 
Finally, for even weaker magnetic fields, the weak turbulence formalism again applies, 
except that now the turbulence is isotropic, i.e. we put fej_ = fen = fe and thus recover 
equation (71). Table [TJ summarizes the essential properties in the three regimes. 

Using up to 2048 3 simulations of decaying MHD turbulence with different initial 
field strength, Lee et al. |187] showed that all three scalings are indeed possible. In 
Figure 12 we show compensated power spectra for three runs with different initial field 
strengths with t>A/u rms w 0.9, 1.3, and 2.0, that are consistent with the regimes of 
Iroshnikov-Kraichnan turbulence, strong turbulence, and weak turbulence, respectively. 



Table 1. Summary of the essential properties of the three regimes of MHD turbulence. 





Iroshnikov-Kraichnan 
(isotropic, sub-equip.) 


strong turbulence 
(critically balanced) 


weak turbulence 
(wave turbulence) 


VA/u lms ~ X 1 


< 1 


~ 1 


> 1 


7"casc 


X~ 2 ta (with fej_ = fey) 


X^ta (= t nl ) 


X~ 2 ^A 


e 


z%k/v A 


4J± 


4 ± k ±/ v Ak\\ 


fe±/fell 


1 


oc k 1 / 3 


— > oo 


E(k x ,k\{) 


(ev A ) 1/2 k- 3 / 2 


e^kl" 3 


(e VA fe||) 1/2 fe: 2 
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Figure 12. Total energy spectra compensated by fc 5 / 3 and averaged over At = 0.5 
(1.5 to 2 turnover times) about the maximum of dissipation for three runs: solid line for 
super-equipartition initial fields (vA/u rms ~ 2.0, dashed for equipartition initial fields 
(va/mi-his ~ 1-3) and dots for sub-equipartition initial fields (vA/u Tlas ss 0.9). The 
three arrows indicate the magnetic Taylor scale. Note that the three spectra follow 
noticeably different spectral laws and possibly different scale-dependence for their time 
scales as well. In all cases the numerical resolution is 2048 3 . Courtesy of E. Lee et al. 
[187] . and Copyright (2010) by The American Physical Society. 



6.5. Dynamo action 

In the absence of an externally imposed magnetic field, it is possible that the field-free 
state is unstable to the dynamo instability, which leads to a conversion of kinetic to 
magnetic energy. If dynamo action occurs, the magnetic field will grow exponentially 
to become dynamically important. The precise outcome regarding energy spectra and 
structure functions is still uncertain, but there is mounting evidence that in the inertial 
range they are similar to those in the purely hydrodynamic case |104[ I196[ I197j ; see 



Figures 13 and 14 However, the largest resolutions obtained in MHD simulations 
so far are still only between 1536 3 [198J and 2048 3 mesh points [187J, and it is not 
necessarily surprising that there is no evidence for a clear bottleneck effect, although 
the spectra have always been seen to be slightly shallower than A; -5 / 3 and closer to 
£r 3/2 [HH QjSJ |200l [2011 HE]. However, it has been argued that, compared with fluid 
turbulence, MHD turbulence is more nonlocal in spectral space |186j . The anticipated 
spectral bump would be more spread out, which might explain the absence of a 
bottleneck at the resolutions available so far, and that much larger resolution would 
be needed to see it. For supersonic MHD turbulence with dynamo action there is 
evidence that the mixed longitudinal structure functions of Politano & Pouquet [28] in 
equation (18) scale linearly with r, provided the Elsasser variables are scaled with a p 1 ^ 3 



factor [2031 [204J. 
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Figure 13. Third order structure function, as denned in equation (16 1, for 
hydromagnetic runs with 512 3 mcshpoints [196] . The inset gives the result for 256 3 
meshpoints. The scaling for transversal structure functions (dotted lines) tends to be 
better than for the longitudinal ones (solid lines). The statistics for the 256 3 runs is 
somewhat better than for the shorter 512 3 runs. Courtesy of Nils E. Haugen 




Figure 14. Structure function exponents for interstellar turbulence simulations (left) 
and fractal dimension of the dissipative structures (right) from simulations. Courtesy 
of Miguel de Avillez [197]. 



In the absence of helicity and with full isotropy, a successful dynamo (positive 
growth rate in the linear regime or finite amplitude in the nonlinear regime) is referred 
to as small-scale dynamo. This refers to the nature of the dynamo process rather than 
just the typical scale of the magnetic field. For example, a small-scale magnetic field 
that is just the result of shredding of an imposed large-scale field is not the result of 
any dynamo process. On the other hand, in the presence of helicity, or with anisotropy 
combined with a mean shear flow, there is the possibility of large-scale dynamo action. 
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Figure 15. Kinetic and magnetic energy spectra in a turbulence simulation without 
net helicity (left) and with net helicity (right) for a magnetic Prandtl number of unity 
and a mesh size is 512 3 meshpoints. Notice the pronounced peak of M(k) at k = k% 
in the case with helicity. Adapted from Refs. [TT] and |205| . respectively. 



6.6. Large-scale and small-scale dynamos 

A typical large-scale dynamo produces magnetic energy on a scale larger than the scale of 
the energy-carrying eddies. A small-scale dynamo is one that generates magnetic energy 
on scales smaller than the scale of the energy-carrying eddies. The difference between 
large-scale and small-scale dynamos is demonstrated in Figure [15] where we compare 
kinetic and magnetic energy spectra of turbulent dynamos with |205] and without [UJ 
helicity. Flows with a large-scale helical pattern of alternating sign, such as the Taylor- 
Green flow |206[ I207j . may be considered as an intermediate case between large-scale 
and small-scale dynamos. 

In the following we use the magnetic Reynolds number, defined analogously to the 



fluid Reynolds number in equation (11 ) by replacing v by r), 

Re M = Wrms/^f- (81) 

The ratio between kinematic viscosity and magnetic diffusivity is referred to as the 
magnetic Prandtl number, Ptm = v/v- The onset of a dynamo is characterized by 
ReM > ReM.critj where Recent is the critical value. An important difference between 
large-scale and small-scale dynamos is the different dependence of ReM.crit on Pr^f. 
Establishing an asymptotic dependence of ReA/ jCrit on Pi M is important because, even 
though the computing power will increase, it will still not be possible to simulate realistic 
values of Ptm in the foreseeable future. Schekochihin et al. [208] have compared the 
results from two independent codes and show that there is a sharp increase of Recent 



with decreasing Ptm', see Figure 16 (where the two quantities are denoted as Rm c and 
Pm). Such a result was first derived analytically |209j . well before it was seen also in 
simulations. 




Figure 16. Dependence of i? cr it on Re. The plot on the left hand side compares the 
results from spectral simulations with those of meshpoint methods. The plot on the 
right hand side shows the results of the most recent simulations with magnetic Prandtl 
numbers down to 0.06 and compares with the case of the Taylor-Green flow |206U207] . 
"JLM" refers to simulations done with the incompressible spectral code written by 
J. L. Maron: runs with Laplacian viscosity, 4th-, 6th-, and 8th-order hyperviscosity 
(resolutions 64 3 to 256 3 ). In this set of simulations, hyperviscous runs were done at the 
same values of 77 as the Laplacian runs, so the difference between the results for these 
runs is nearly imperceptible. "PENCIL" refers to weakly compressible simulations 
done with the Pencil Code: runs with Laplacian viscosity, 6th-order hyperviscosity, 
and Smagorinsky large-eddy viscosity (resolutions 64 3 to 512 3 ). Courtesy A. A. 
Schekochihin [2135] , 



The reason for the increase of Recent with increasing Re has been explained by 
Boldyrev and Cattaneo |210] as being related to the fact that when Re^ < Re, the 
resistive scale (i.e. where the magnetic power spectrum peaks in the kinematic regime) 
shifts from the dissipative subrange into the inertial range. In the inertial range the 
velocity field is no longer smooth, but it is rough in the sense that the exponent £1 (see 



2.4) in the scaling of velocity differences over distance r is less than 1 |210j . For £1 < 1, 



the velocity field becomes non-differentiable in the sense that velocity gradients diverge 
like r'* 1-1 . The smaller (1, the rougher the velocity field, while (1 = 1 corresponds to a 
smooth velocity field. 

More recent work |21R 1212] suggests that the threshold for small-scale dynamos 
is particularly high only in the range 0.06 ^ Prjv/ ^ 0.2, because then the resistive 
scale lies within the range where the kinetic energy spectrum shows the bottleneck 
with (1 — > 0, corresponding to an extremely rough velocity field with very large critical 
magnetic Reynolds number. However, when Pr A / ^ 0.06, the resistive scale lies beyond 
the bump of the bottleneck, i.e. well inside the inertial range, and there the critical 
magnetic Reynolds number is again somewhat smaller. Resolving this issue conclusively 
requires a numerical resolution well in excess of 1024 3 meshpoints, as well as long run 



Astrophysical turbulence modeling 



42 



10.00 



l.oo 



0,10 



Dl 



10.00 



1,00 ■ 



0.10 - 



O.Ot L 



1 



1 'd 



Pr w =O.0l 



]00 





1 


■ — 






















-_ 

■-. 


































: ✓ 






\ \ 




: / 
.' 






\ ^ 


■ 


( 











13 



ino 



10.00 




V \ T o.io 





1 






— R ■ : 

; 




1 












^" -"■ 












■ •._ •-. 






(/ 




''■ ■ - ■. 










X ^ 






















\ \ 1 


■ /■ 








A \ : 

\ n : 










\ v 













10.00 



1.00 



0.10 



0.01 



I 



Pr^O.001 



100 



' X. 

\ */ 


■ ■ | ■ ■ ■ — ■ — ■ — ■ ■ ■ i ■ — : 

" v N E{k)^i^ \ 
\ - - 


: ' k 
■ f K Ji 

r 


\ \ 

\ ^ 


10 100 



Figure 17. Compensated kinetic and magnetic energy spectra in the saturated 
regime for Pr M = 1 with Re = 450, Pr M = lO" 1 with Re = 1200, Pr M = 10~ 2 
with Re = 2300, and Ptm — 10" 3 with Re = 4400. The spectra are compensated 
by ey 2 / 3 fc 5 / 3 , where e? is the sum of kinetic and magnetic energy dissipation rates. 
The ohmic dissipation wavenumber, k„ = (e^//?? 3 ) 1 ^ 4 , is indicated by an arrow. The 
viscous dissipation wavenumbers are 180, 290, 350, and 430 for Prjv/ = 1, 10 _1 , 10 -2 , 
and 10 -3 , respectively. Adapted from Ref. 205 . 



times, which is only now beginning to become feasible. We may therefore expect further 
developments in this area in the near future. 

If there is large-scale dynamo action, the magnetic field grows preferentially at 
scales large compared with the energy-carrying scale. This process is non-local in 
spectral space [213] , although it has also been shown that an externally applied magnetic 
field produces mainly local interactions |214j . On the other hand, large-scale dynamo 
action depends on velocity and magnetic field correlations at the energy-carrying scale 
(rather than the resistive scale). The onset of this type of large-scale dynamo action is 
essentially independent of Pr M and occurs when Re^ > Recent ~ 1- The independence 
of the saturation strength of the large-scale dynamo on the microscopic resistivity is 



demonstrated in Figure 17, where we show spectra of kinetic and magnetic energies for 
different values of PrM- 
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6. 7. Turbulent convection and stratification 

In certain layers of a star the opacity of the gas can become so large that the 
energy flux can no longer be transported by radiative diffusion, but by convection. 
A phenomenological theory called mixing length theory allows one to make reasonable 



estimates for the expected turbulent velocity. As mentioned in § |3.2[ the convective 
energy flux is approximately equal to the pu^ ms . This gives a good estimate for the 
convective velocity in a star. 



The basics of the convection instability was discussed in § A necessary condition 
for convection is that the specific entropy decreases with height, i.e. iVg V < 0; see 



equation (47). In addition, viscosity and thermal diffusion have to be small enough 
compared with the height of the unstable layer, d, and the Brunt-Vaisala frequency, 
iV B v- This is quantified by the Rayleigh number, 

Ra=-(-A^ v ) (82) 
u X \ Jo 

which has to be above a certain critical value for the onset of convection. Here, 
the subscript refers to the requirement that the specific entropy gradient has to 
be calculated for the associated hydrostatic equilibrium solution, and not for the 
already convectively unstable solution. Such solutions are not normally presented in the 
literature. Also, the thickness of the outer layers of the Sun would be much smaller in 
the hydrostatic reference state. It is therefore not common to quote Rayleigh numbers in 
astrophysics, except in idealized simulations whose hydrostatic reference solutions tend 
to be polytropes where the initial density is related to the initial temperature via p ~ T n , 
where n is the polytropic index. Unlike the incompressible case, where the Rayleigh 
number is based on the background gradient of temperature, in the compressible case it 
is based on the gradient of specific entropy for the associated hydrostatic solution [52J . 

If the value of the Rayleigh number is increased sufficiently beyond the critical value, 
the flow becomes turbulent. Simulations of turbulent convection have been provided by 
many different groups, both in the incompressible approximation |215[ I216[ 1217] as 
well as in the fully compressible case |148[ 1218] 1219] . Typical Rayleigh numbers that are 
currently reached in simulations are around 10 6 . With rotation the onset of convection is 
delayed correspondingly, which enables one to reach somewhat larger Rayleigh numbers 
in such cases. 

The Nusselt number is another commonly used quantity in incompressible and 
laboratory convection. In that case it gives the ratio of the total heat flux to that 
transported by heat conduction alone, using the same boundary conditions. However, 
unlike laboratory convection, where the temperatures at top and bottom are usually 
kept fixed, in many compressible simulations with a polytropic background solution 
the energy flux at the bottom is actually prescribed. One compares therefore normally 
with the radiative solution with a linear temperature profile that has the same top and 
bottom temperatures as the convective solution. One also subtracts out the flux that 
is transported by the adiabatic stratification alone |151j . Again, this value is nowadays 
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Figure 18. Snapshots of B y in the early phase (left) and saturated phase (right) 
of the dynamo. The sides of the box show the periphery of the domain whereas the 
top and bottom slices show B y at top and bottom of the convectively instable layer, 
respectively. Courtesy Petri J. Kapyla [226 . 



not normally quoted for compressible simulation. For many purposes, a more useful 
characterization of the turbulence is the resulting value of the Reynolds number. 

Another important difference to laboratory convection is the absence of boundaries 
in astrophysical convection. Convectively unstable layers are the result of a particular 
dependence of opacity on temperature and density. This has frequently been modelled 
by using prescribed spatial profiles of the radiative conductivity. In this way one can 
model convection in an unstable layer, sandwiched between two stable layers [220] . This 
makes the dynamics near the transition layer softer and allows the flow to overshoot 
into the stably stratified layers. This leads to the excitation of gravity waves in the 
stably stratified layers [2201 [22H [2221 12231 [221 1225] - 

Convective flows can well support dynamo action. As an example we mention here 
the result of a convection simulation with horizontal shear which leads gradually to the 
development of a large-scale magnetic field [226] . A result of such calculations is shown 
in Figure 18, where we visualize the toroidal field component at an early time when only 
small-scale fields have been produced, and at a later time when also a large-scale field 
is present. 

The presence of large-scale fields is often characterized by energy spectra. However, 
because of stratification it only makes sense to look at horizontal spectra taken at a 
specific depth. If the mean magnetic field depends mainly on depth, the horizontal 
magnetic energy spectra will peak at wavenumber zero, which can only be seen if 
one plots the spectral energy versus linear wavenumber; see, for example, Fig. 12 of 
Ref. 
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6.8. Global hydromagnetic dynamo simulations 

Simulations of global convection have demonstrated the generation of differential 
rotation and magnetic fields [22 7\ I228j . However, with parameters relevant to the Sun 
such models have not yet produced large-scale magnetic fields similar to those in the 
Sun [551 1121 j . This is plausibly explained by the relevant dynamo numbers for coherent 
large-scale field generation being still too small. In that case, only small-scale magnetic 
fields are generated, while the threshold for large-scale field generation has still not been 
reached. This is different when the rotation rate of the sphere is increased to several 
times the solar value [229J. As an example we show here the results for a sphere that 
has a stratification similar to that of the Sun, but it is rotating about 3 times as fast; 
see Figure [T9| 

The rapid rotation is primarily responsible for producing the typical convection 
patterns that are elongated in the direction of the rotation axis. This effect is especially 
obvious at low latitudes, outside the inner tangent cylinder, i.e. the cylinder that is 
tangent with the bottom of the convection zone. The resulting convection pattern is 
often referred to as banana cells, a concept that was widely discussed in the 1980s [230] , 
but there has never been observational evidence supporting this type of flow pattern for 
the Sun. Banana cells occur as a consequence of rapid rotation, which is also responsible 
for cylindrical angular velocity contours. Although this does not apply to the Sun, it 
may well apply to some stars that rotate much more rapidly than the Sun. 

Simulations of rapidly rotating convection [229] l5*o] show that in the region with 
strong banana cell convection, there is strong large-scale dynamo action with pronounced 



toroidal flux belts on one or both hemispheres; see Figure 19 This is partially 
reminiscent of the magnetic activity in the Sun, although it would be premature to 
draw any conclusions from this given that at present there is no explicit evidence of 
banana cell convection in the Sun. 

We mention here another line of research. Instead of convection driving the flow 
one can apply an artificial forcing function. This has the advantage of producing a flow 
pattern whose typical scale can be controlled. In particular, it is possible to achieve 
turbulent scales that are small compared with the radial extent of the domain, so as to 
produce a well-defined scale separation [2311 1232] . With such simulations it has been 
possible to to focus entirely on the nature of the dynamo in spherical shell geometries 
and to isolate its physics from many other effects that may still be important. It 
turns out that even in the absence of global shear, oscillatory large-scale fields can be 
generated |232[ 1233] . Such solutions show equatorward migration and are quite different 
in nature from oscillatory solutions of aVt type. It is well possible that these solutions 
have nothing to do with those in the solar dynamo, but it serves as a reminder that the 
variety of possibilities can be much larger than what is usually discussed. 
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Figure 19. Toroidal and radial magnetic field (first and second row) together with 
radial velocity (bottom row) near the top of the convective shell (left column at 
r/R = 0.95) and in the middle (right column at r/R = 0.85). The mean magnetic field 
is approximately antisymmetric about the equator. The radial velocity shows flows 
patterns elongated along the rotation axis (so-called banana cells). The resolution 
is 96 x 256 x 512 mesh points or collocation points in the radial, latitudinal, and 
longitudinal directions, respectively. The magnetic Reynolds number based on the 
thickness of the convective shell and without dividing by 2ir is 86, and the Coriolis 
number, i.e. the ratio of vorticity from the mean rotation to the rms vorticity of the 
turbulence, is about 3. Courtesy Benjamin P. Brown [229 . 



6.9. Interaction between convection and shear 

Simulations of rotating convection in spherical shells demonstrate that there is 
equatorward acceleration of the mean flow. This phenomenon is generally referred to as 



differential rotation and will be discussed in more detail in § |7.3[ In addition, one sees 
that the convection pattern itself moves differentially across the surface. However, a 
more detailed inspection reveals that at the equator the convection pattern can actually 
move still somewhat faster than the mean flow. This has been revealed both by linear 
theory [234, 235l l236j and by nonlinear simulations (237], and may explain a phenomenon 
seen at the solar surface which shows that magnetic tracers move at speeds faster than 
the speed of the plasma. In fact, even very young sunspots tend to move not only faster 
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than the plasma at the surface, but they move also faster than the gas at any other 
place in the Sun, as seen by global helioseismology; see Figure 4 of Ref . [238J . 

There is at present no universally accepted theory for the enhanced rotation speed 
of magnetic tracers on the Sun. It has however been pointed out that the enhanced 
pattern speed of magnetic tracers might be understandable if the observed magnetic 
field (including that responsible for producing sunspots) was generated in a layer not 
too far below the surface [239]. This proposal would be in conflict with the generally 
adopted view according to which the magnetic field responsible for the solar cycle is 
generated near or even below the bottom of the convection zone of the Sun. 

6.10. Granulation, convection, and solar abundances 

Simulations of solar granulation have reached a high level of realism and have 
proved to be a viable and feasible alternative to earlier one-dimensional models for 
calculating diagnostic spectra in visible light. Strictly one-dimensional models always 
needed to incorporate ill-determined parameterizations of what is known as micro 
and macro turbulence. New realistic three-dimensional simulations of solar convection 
021 EH EH GUI EIH [2321 EE] lead to diagnostic spectra that can be fitted to observed 
spectra without invoking these ill-known parameterizations. The use of 3-D models also 
results in abundances derived from different spectral features (e.g. molecular and atomic 
lines) being more consistent. 

Initial efforts to derive updated solar abundances based on 3-D models resulted 
in new abundance estimates for the heavier elements in the Sun that were as low as 
only 60% of previous estimates |244j . It should be noted, however, that even though 
these abundances are often referred to as "3-D abundances" , 3-D effects were not the 
main cause of the systematic lowering of the abundance estimates, which were instead 
a combined result of updated oscillator strengths, different line fitting procedures, and 
choices made when estimating collision cross sections important for non-LTE corrections 
for some spectral lines. This was elucidated by an independent analysis by a different 
group [501 1243] , who confirmed that 3-D effects improve the consistency but do not give 
rise to a significant systematic abundance effect for the important heavy elements. 

The abundances of the heavier elements determine the opacity of the gas and 
thereby the detailed radial structure of the Sun. On the other hand, the radial 
dependence of the sound speed and density in the Sun can be determined independently 
through helioseismology [591 EOl EH I245| 1246] . and helioseismology can thus provide 
important constraints on the heavy element abundances in the solar interior. (It may in 
the future be possible to also determine the Sun's deep interior composition by exploiting 
neutrinos from the CN cycle and the p-p chain to determine the primordial solar core 
abundances of C and N at an interesting level of precision |247j .) In the convection 
zone the gradual ionization of carbon, nitrogen, and oxygen with depth influences the 
equation of state, and helioseismic measurements of the effective ratio of specific heats 
of the gas can thus provides constraints on the abundance of these elements also there 
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[6211218]. 

The significant downward revision of solar abundances proposed in Ref. [242] and 
even the somewhat more moderate revisions proposed more recently by the same 
group [51] turned out to be difficult to reconcile with observational constraints from 
helioseismology, despite many different attempts to do so; cf. Ref. [62J and references 
therein. However, the downward revisions recommended by [50] are only about half 
as large and are in fact consistent with helioseismic estimates of the heavy element 
abundance in the solar convection zone, Z = 0.167; see Table 2 of Ref. [248J. 

Due to gravitational settling the abundances of all elements differ somewhat 
between the convection zone and the radiative interior [249J. Because of rapid mixing 
the abundance levels are constant in the convection zone, but below the convection 
zone the chemical abundances vary with radius in a manner that is influenced by how 
turbulence in the convection zone generates weak overshooting motions in the radiative 
zone, which result in a slow mixing over depth of chemical elements [250J. 

There was always a small departure in sound speed between models and helioseismic 
observations in a narrow region just below the convection zone. With the revised 
abundance estimates by [242] this departure increased from about 0.3% to about 1.2% 
[62] 1251] , while with the abundances recommended by [50] the discrepancy is of the order 
0.6%. Even the smallest of these discrepancies is many times larger than the helioseismic 
measurement uncertainties, and one should thus worry less about the particular size of 
the discrepancy in any one case, and more about the very existence of the discrepancy. 
In general terms, the lack of a detailed quantitative understanding of the overshoot of 
convection below the bottom of the convection zone and the associated slow mixing 
seems to be a likely reason for the discrepancy [BT] . 

An important additional observational constraint on slow mixing below the 
convection zone comes from the depletion of Lithium in the Sun. Lithium is destroyed 
at temperatures that are reached about one pressure scale height (corresponding to 
about 1% of the solar mass) below the convection zone, and the observed depletion (a 
factor of about 160) implies that mixing down to that temperature takes place on a 
time scale considerably shorter than the age of the Sun, but still very large compared to 
convection zone turnover times [BTJ 1252] . Lithium depletion in other stars is now known 
to be essentially consistent with the behavior expected from the differences in age and 
structure deduced from standard stellar evolution theory [252J. 

6.11. Turbulence from the magneto-rotational instability 

In the presence of shear and rotation, the slow magnetosonic waves develop a long 
wavelength instability, where u 2 < for v\k 2 < 2qQ 2 . Here, q = — dlnf2/dlnr 
quantifies the radial gradient of the angular velocity. This is called the magneto- 
rotational instability (MRI). It is particularly simple to analyze if the magnetic field 
is vertical, in which case the instability is purely axisymmetric. However, the same 
growth rates are obtained in the nonaxisymmetric case, if B points in the streamwise 




Figure 20. Toroidal magnetic field component displayed on the periphery of the 
computational domain (color coded). The size of the box is 2ir in all three directions 
and the mesh size is 512 3 meshpoints. The gas is isothermal with a constant sound 
speed of c s = 5f2/fci. Viscosity and magnetic diffusivity are v = 77 = 2.5 x \{)^ i Q,/k\. 
Note that the magnetic field develops long thin structures aligned at some angle relative 
to the toroidal direction. Adapted from Ref. [253] . 

direction [86] . 

In the axisymmetric case the instability takes the form of so-called channel flows. 
In three dimensions the flow experiences strong shearing and hence small length scales 
in the cross-stream direction. This leads to the flow breaking up into what we loosely 
call fully developed turbulence. An example of such a flow is shown in Figure [20| 
where periodic boundary conditions have been used in the vertical and azimuthal 
directions, and shearing-periodic boundary conditions in the cross-stream direction. No 
net magnetic flux has been applied |253] . Numerical simulations show, however, that the 
MRI is no longer excited when the magnetic Prandtl number is less than a critical value 
of the order of unity [2541 1255] . A strong sensitivity on the magnetic Prandtl number has 
also been found for magnetized Taylor-Couette flow [256] 12571 1258] 1259] 1260] . At present 
it is still unclear whether there is a real problem with the MRI in accretion discs when the 
magnetic Prandtl number is small. This issue may well be connected with the difficulty 
to excite small-scale dynamos at low magnetic Prandtl numbers [208, 209j 12101 1211] 1212] . 
On the other hand, astrophysical dynamos are large-scale dynamos, and they do not 
suffer from that particular difficulty [205] 1261] . It would therefore be important to 
perform new MRI simulations in cases where large-scale dynamos are possible, i.e. in the 
presence of vertical density stratification, which can then lead to an a effect [2621 |263j . 
In another recent study it has been shown that even without stratification, large-scale 
dynamo action is possible when pseudo- vacuum boundary conditions are used at top and 
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Figure 21. Visualization of the logarithmic density of an accretion torus around a 
black hole. Courtesy of John F. Hawley |277j . 



bottom of the rotating shearing box [264J. A similar generation of mean fields has also 
been found without rotation |239l 1265} 1266] I267j . Possible candidates for explaining the 
origin of large-scale fields in this case include the incoherent a-shear dynamo |268[ 269J 
and the shear-current effect |270] [271J. For the latter effect to work, it is necessary that 
one of the off-diagonal components of the magnetic diffusivity tensor has a suitable sign, 
which may however not be the case [2721 12671 12731 1274] . 

The MRI is generally thought to be responsible for driving turbulence in accretion 
discs, where q = 3/2. A more accurate representation of accretion discs is obtained 
with the inclusion of vertical and radial density stratification. The former case can 
be treated within the shearing-box approximation |275l 1276] while the latter requires 
a global treatment |277[ 1278} 1279} 1280] . In Figure 21 we show a visualization of the 
logarithmic density of an accretion torus around a black hole 

An important diagnostic quantity of accretion disc simulations is the dimensionless 
turbulent disc viscosity, «ss — v t /c s H. Here, the subscripts SS refer to Shakura and 
Sunyaev |281j . who employed this parameterization of turbulent viscosity z/ t in terms 
of local sound speed c s and pressure scale height H. In the simulations, z/ t is normally 
estimated by the mean total horizontal stress, IIr^ = b^b^j p — puru^, divided by the 
mean rate of strain resulting from the differential rotation, pRdQ / dR, where cylindrical 
coordinates, (R,(p,z), have been employed. 

In comparison with local shearing box simulations, an important difference is that 
global simulations are capable of producing about 10 times larger values of «ss- This 
is an immediate consequence of the larger field strength in global simulations rather 
than a difference in the intrinsic properties of local versus global disc simulations |277] . 
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Figure 22. Simulations showing the total column density (gas plus all particle sizes) 
in the horizontal plane. The insets show the column density in logarithmic scale 
centered around the most massive cluster in the simulation. Time is given in terms of 
the number of orbits after turning on self-gravity. Courtesy of Anders Johansen. 



Another important outcome of global disc simulations is the fact that Ur^ is finite at 
the innermost marginally stable orbit. This is a property that is not normally taken 
into account in analytic models and continues to be debated in the literature [2821 1283] . 

A number of new simulations have emerged in recent years. A major step was 
the combination of dust dynamics with self-gravity in the shearing box approximation 
[284} I285j . One of the remarkable results is a rapid formation of nearly Earth-sized 



bodies from boulders (Figure 22). Even though the mass of what one might call 
protoplanet is growing, this body is also shedding mass during encounters with ambient 
material as it flows by. One might speculate that what is missing is the effect of radiative 
cooling of the protoplanet. This would allow the newly accreted material to lose entropy, 
become denser, and hence fall deeper into its potential well. 

The main reason the simulations presented in Ref. |284j produce rapid growth 
is connected with the occurrence of sufficiently strong compressions caused by the 
turbulence. Once the compression is strong enough, self-gravity takes over and leads to 
a fully developed nonlinear collapse. 

6.12. Effects of thermal and gravitational instabilities 

A thermal instability may arise if a cooling term, A(T), and a heating term, T(T), are 
included on the right hand side of the energy or entropy equation, i.e. 

pT^ = ... + pT(T)-p 2 A(T). (83) 




Figure 23. Visualization of InT on the periphery of the box at different times, for 
v = x = 5 x 10~ 4 Gyr km 2 s~ 2 and 256 3 mesh points, (p) ps 1.74 x 10~ 24 gcm~ 3 , 
(p) 24.2 x lCT 14 dyn, and (T) ps 8200 K. Here Q = 100 Gyr" 1 and S = -fl. For this 
run k-p/ki — 32 and n p /(c s kp) = 1.5. The growth rate is about 190 Gyr^ 1 , which is 
somewhat larger than for the corresponding non-shearing run. Note that the initially 
produced structures are quickly sheared out. Adapted from |163j 



It is convenient to abbreviate the combination of the two terms on the right hand side 
by pC, where £ = T — pA. This allows us to state a sufficient condition for stability 
[To31l2gE] 

{w) >0 ( stabilit y)- ( 84 ) 

This means that when the temperature is increased, the corresponding cooling increases, 
bringing the temperature down again to the original value. In the presence of thermal 
diffusion, with = —K'VT ^ 0, the system can always be stabilized at small scales, 
i.e. for large wavenumbers, where eventually the thermal diffusion rate becomes faster 
than the cooling rate. For T = const and A oc T@, the dispersion relation u(k), is on 
sufficiently large scales (small wavenumbers) of the form |163l 1286] . 



io = c iso k^l - (3~\ (85) 

where Ci SO = c s / y/j is the isothermal sound speed. Evidently, for j5 < 1 sound waves 
become destabilized (u becomes imaginary). 

Numerical simulations |162l 1163] have not been able to confirm alternative findings 
[287] that the thermal instability can lead to sustained turbulence. This is demonstrated 
in Figure 23, which shows (here in the presence of shear) that the thermal instability 
leads to the development of patches with low temperature (100 K compared to 10,000 K 
outside those patches), but over time these patches merge until eventually a stable 
equilibrium is reached where a few big patches continue to coexist. 

Another instability where sound waves are destabilized is the Jeans instability. Here 
the dispersion relation can be written in the form |288j [289, 290, 1291] 

to 2 = c 2 k 2 - AnGp, (86) 

where p is the local density of the gas. So, again, large scales become unstable. In an 
asymptotically thin layer such as an accretion discs or galaxies, the dispersion relation 
becomes [292J 

uj 2 = c 2 s k 2 - 47rGS|fc|, (87) 
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where E is the local surface density. In the context of local accretion disc models, 
simulations suggest that this process can indeed lead to sustained turbulence [871 293, 
12941 1295] . In simulations of star formation (TEH [2961 EHTJ ESS], the Jeans instability 
leads to a continuous production of gravitationally bound structures corresponding 
to protostars. The stars that form have a broad distribution of masses, determined 
mainly by the statistics of mass fragmentation in supersonic MHD-turbulence |166j . 
The fraction of stars that are heavier than about 8 solar masses eventually (after a 
delay of up to a few tens of Myr) explode as supernovae. These supernovae contribute 
to sustaining the turbulence in the interstellar medium that ultimately causes additional 
generations of stars to be born [721 If 97] . 



6.13. Supernova- driven turbulence 

Interstellar turbulence is an example of astrophysical turbulent flows where the driving 



is usually modeled by a distributed body force. As discussed in § |3.4[ the blast 
waves of supernova explosions provide energy input to the surrounding gas. These 
explosions drive gas flows with temperatures of around 10 s K, but they also lead to 
strong compressions where the gas cools rapidly to about 10 4 K. When the temperature 
is between 100 K and 10 4 K the gas may, depending as details of the cooling curve A(T), 
be thermally unstable [286]. This contributes to keeping the gas in the interstellar 
medium preferentially in one of two distinct temperature regimes (the so-called cold 



and warm phases; see §6.12). The hot phase at temperatures > 10 K is a direct 
result of heating by supernova explosions combined with a low cooling efficiency of the 
interstellar medium at that temperature. This is also borne out by various simulations 
[2991 [3001 BOD • Simulations show that the filling factor of the hot gas (T > 10 6 K) grows 
with height from 0.2-0.3 at the midplane to about 0.5-0.6 at a height of about 300 pc 
[299J . However, this result depends on the degree of correlation of supernovae in space 
and can reach 0.6 at the midplane for completely uncorrelated supernovae, as in an early 
analytic model [302] . Simulations have also been able to demonstrate that significant 
amounts of vorticity are being produced if the flow is sufficiently supersonic and if 
the baroclinic term is important [303, 1304] . The presence of vorticity is advantageous 
for dynamo action; in fact, no dynamos have yet been found when the turbulence is 
irrotational 305 . 



There is now mounting evidence that for large Mach numbers the energy ratio 
of compressive to solenoidal velocities approaches 1/2 [1661 H791 13061 13071 13081 [309J. 
This can be explained if the mean square values of longitudinal and transversal velocity 
derivatives were equal, i.e. (u 2 x ) = (u xy ). Assuming isotropy and that mixed terms 
cancel, this implies ((V ■ u) 2 ) « 3(ul^) and (uj 2 ) ~ G(u 2 y ), giving a ratio of 1/2 [307] . 
Whether or not this behavior is really universal needs to be seen. In the papers listed 
above the turbulence was forced with a substantial solenoidal component, so the issue of 
vorticity production was not addressed. In the following we discuss the opposite limit, 
where only compressive modes are driven and where no vorticity is produced. 




Figure 24. Visualization of lap on the periphery of the box at different times, for 
kiR = 0.2, Rc = 50, and 512 3 mesh points. Note that in the fully developed state 
individual expansion waves can hardly be recognized. Adapted from Rcf. [305J . 



6.14- Irrotational turbulence 

Turbulence is usually thought of as being an ensemble of interacting eddies. If one 
associates eddies with vortices, then "irrotational" turbulence must be a contradiction 
in terms. Nevertheless, irrotational turbulence can be regarded as an idealization that 
can serve its purpose in illustrating the difference to regular (vortical) turbulence. 



Irrotational turbulence means that w = V x « = 0. As explained in § the 
u x oj nonlinearity is absent and the only nonlinearity comes from the \u l term. This 
causes a significant modification of the turbulent cascade, which is one of the reason 
why irrotational turbulence may be a contradiction in terms. Because of compressibility, 
however, vorticity can in principle be generated via the viscous term. Taking the curl 



of • r in equation (33), and assuming v = const, gives 




V x jJ-V • rj = vV 2 u> + V x [2vS ■ V ln(pv)] . (88) 
Here, the first term vanishes if oj = 0, but the second term does not. As mentioned 



m 



6.13, simulations show that this term remains small in the limit v — > [305J. 



In Figure 24 we show visualizations of the logarithmic density in a simulation, which 
shows that the initially highly ordered expansion waves turn rapidly into a complicated 
pattern. The flow is here driven by a forcing function / = — V0, where is a scalar 
function consisting of randomly placed Gaussians that change in regular time intervals, 
At, such that Atu rras k{ ~ 0.25. 

Given that viscosity always perturbs the zero vorticity state slightly, and because 
the vorticity equation is analogous to the induction equation, one must ask whether a 
small initial vorticity could increase owing to an instability. However, at the Reynolds 
numbers achieved so far in simulations, neither vorticity nor magnetic fields have 
been found to increase spontaneously |305] . The suggestion that purely irrotational 
turbulence cannot produce dynamo action may be related to the finding that in vortical 
supersonic turbulence the critical magnetic Reynolds number for small-scale dynamo 
action shows a "bimodal" behavior with Mach number: for Mach numbers below unity 
the critical magnetic Reynolds number is about 35 to 40, and above unity it is about 70 
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to 80 |307] . Note, however, that the flow is here not purely irrotational, and that the 
ratio of ((V ■ u) 2 ) and (uj 2 ) is about 1/2; see the discussion in §6.13 



The results concerning vorticity production may be of relevance for other flows that 
can be described by spherical expansion waves. One example concerns phase transition 
bubbles that are believed to be generated in connection with the electroweak phase 
transition in the early universe [3101 1311] . Here the equation of state is that of a 
relativistic fluid, p = pc 2 /3, where c is the speed of light. Thus, again, there is no 
baroclinic term and no obvious source of vorticity. However, the relativistic equation of 
state may be modified at small length scales, but it is not clear that this can facilitate 
significant vorticity production. 



7. Collective effects of turbulence 



In this section we denote the velocity by a capital U. Overbars indicate averages over 
one or two coordinate directions. They are not therefore regarded as spatial filters that 
are often used in the theory of large-eddy simulations [3121 1313] . The definition of 
averages in terms of coordinate averages is convenient for interpreting simulation data. 
Other definitions of averages are possible. In analytic studies ensemble averages are 
commonly used. Departures from these averages are denoted by lower case symbols, 
i.e. u = U — U and b = B — B denote the fluctuating components of the velocity and 
magnetic field vectors. We discuss the properties of various correlators such as UiU], 
Uibj, and b, t bj. 

In general turbulence is non-isotropic. This can lead to the possibility of non- 
trivial components of the correlations tensors ujllj, u^bj, and bibj. The effects of these 
correlations on the evolution of the mean flow, U, and the mean magnetic field, B, or the 
mean passive scalar concentration, C, are referred to as collective or mean-field effects 
of the turbulence. Even in the special case where UjUTj is a diagonal tensor there is at 
least the phenomenon of turbulent diffusion, which will now be illustrated in connection 
with the passive scalar field. 



7.1. Turbulent passive scalar diffusion 

The relevant dynamics comes from the nonlinearity. In order to keep the discussion 
simple, we neglect the diffusion term. The evolution equation of the passive scalar 
density per unit volume, C = pO, is then 
DC d 



dt dxj 



(CUj) , (89) 



cf. equation (43). Again, we define the average concentration per unit volume as C and 



write C = C + c. The evolution equation of C is obtained by averaging equation (89) 
i.e. 

dC d 



dt dxj 



CUj + cuj) . (90) 
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The problematic term here is cuj, and the hope is that it can be expressed in terms of 
mean fields such as C and U. 

In order to derive an expression for cu], we consider its evolution equation, 

8 

—c Uj = c Uj + cuj, (91) 

where dots denote partial time derivatives. The evolution equation for c is obtained by 
subtracting equation (90) from equation (89), which yields 

i=-4^ + ^ +J ^)- (92) 

where i\T- = cuj — cu] denotes nonlinear terms. In the absence of rotation, shear, 
viscosity, or other linear effects, the momentum equation takes the form 

Nf. (93) 

Assuming incompressibility and no mean flow, U = 0, we have 



d _}h = at(«) 
dt 3 



9 ~ ( c ) dC 

dt CUj ~ H " dx ' ' 



where k-j = UjU] and r( cu ) = — [V-ZV^jit + cNW denotes a triple correlation term. 

Clearly, in the statistically steady state the two terms on the right hand side of 
equation (94) must balance to zero, suggesting that cannot be neglected, as is 

assumed in the commonly used first order smoothing approximation, when it is applied 
to the case of vanishing diffusivity (kq = 0); see Ref. |314] for a more detailed discussion. 
When Kg is large, the microscopic diffusion term involving KqS/ 2 9 in equation (43) or 
KcV 2 C in equation (89) needs to be restored. Since it is applied to the small-scale 
field with typical wavenumber k{, the inclusion of the Kg term corresponds essentially 
to adding —KekfcUj on the right-hand side of equation (94). (This can be treated more 
accurately in Fourier space; see Ref. |315j for a corresponding treatment in the magnetic 
case.) 

The closure assumption used in the r approximation consists of the assumption 
that the triple correlations can be expressed in terms of the quadratic correlation, i.e. 

Tj = —cuj/t (closure assumption). (95) 



Inserting this into equation (94) yields 



d \ (a dC 

3 




where k\^ = tk\j corresponds to the usual turbulent diffusivity. This equation shows 
that, in the statistically steady state, there is a flux of passive scalar concentration in 
the direction of the negative gradient of C. Note that the effect described here works 
also when the turbulence is isotropic, i.e. when UjU] = \8ijU 2 . In that case we have 
Ky = K^Sij, where is the scalar turbulent diffusivity of the mean passive scalar 
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concentration. By assuming r = (urms&f)" 1 we obtain « t = lurms^f 1 - On the other 
hand, if Pe <C 1, k[ is small and increases linearly with Pe such that n t = ^Peu^gk^ 1 . 

The effect discussed above is known as turbulent diffusion. It is a very basic effect 
that characterizes an enhanced diffusion experienced by the mean concentration. It is 
present whenever the typical scale of the mean field is large compared with the scale of 
the turbulence. This is the requirement of scale separation that needs to be made in 
order for a multiplicative relation in terms of the product of n t and VC to be valid. 
On the other hand, if the scale of the turbulence is comparable with the system size, 
a local connection between flux and gradient becomes invalid, and nonlocal expressions 
must be considered [316J. 

Let us now contrast the r approximation with the first order smoothing 
approximation, where equation (92) is still used, but the Nj term is now neglected. 
Again, assuming U = and integrating in time, we have 

l{x 
dx 



c(x, t) = - f d °l X ; tn> Uj{x, t') dt'. (97) 







Thus, 



dC(t') 



cui = - / Ui(t)uj(t') — dt', (9£ 



where we have dropped the common x dependence of all variables for clarity. This 
expression would be identical to equation (96) in the special case where 



Ui(t)uj{V) = -?ZpJexp[-(i - t')/r] for t > t'. (99) 

This type of agreement is restricted to the simplest case when there is no contribution 
from the momentum equation. Examples where such agreement is lost include cases 
with rotation or shear, as well as analogous cases with magnetic field where there can 
be contributions from the Lorentz force [TH 1317] . 

The concept of turbulent diffusion carries over to vector fields such as the velocity 
itself and the magnetic field. In these cases one talks about turbulent viscosity, z/ t , and 
turbulent magnetic diffusivity, r} t . The relevant correlations are then UjWj and Uibj that 
are being expressed in terms of negative gradient terms, i.e. 

fdUi dUA , , 

^=-"'U; + *t)' (100) 

dB 

^=-m^- (101) 



This last formula is quite analogous to the passive scalar case discussed in equation (96), 

^=-4 c) |^, (102) 

where we have dropped the time derivative of UjC. The term on the right-hand side of 
equation (100) is similar to the expression for microscopic diffusion, see equation (32). 
The correlation that enters in the mean induction equation is 



Si = (u x b)i = e ijk Ujb k = -r/t(V x B)t = -r] t /j, Ji, (103) 
which gives a contribution similar the microscopic diffusion term in equation (58). 
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7.2. The a effect 

Turbulence does not always act just diffusively. There can be non-diffusive effects, 
especially if the turbulence lacks local isotropy or at least parity invariance. If the flow 
is statistically non-mirrorsymmetric (for example helical) interesting effects can occur 
in connection with the evolution of the mean magnetic field. In particular, there are 
terms proportional to the mean magnetic field itself, i.e. 

zTy<b = aB- r] t V x B. (104) 

This is the famous a effect pj)J I318[ 1319] that is often invoked in order to understand 
the generation of large-scale magnetic fields in astrophysical bodies. The possibility 



of magnetic field generation can be seen by inserting equation (104) into the mean 
induction equation, 
dB 



V x (u x b - rjV x B) . (105) 

One can look for solutions proportional to exp(ifc • x + At) and finds that 

A = ±ak — (t) + i] t )k 2 , (106) 

where k = |fe|; see, e.g., Ref. [H] for details. This shows that exponentially growing 
solutions exist on sufficiently large scales, i.e. on sufficiently small wavenumbers, 
k < a/r]T- Here we have introduced the total magnetic diffusivity, t/t = f] + T) t . 

Although this topic has reached text book level already several decades ago 
[T5] I318[ 1319] . it continues to be a field of intense research - especially with regards 
to nonlinear feedback. Basic aspects of the a effect can best be explained in the context 
of isotropic turbulence. In that case the following expression for a has been derived 
[320 | I32T | 1322] 

a = -±t(«.v> + ±T<j •&>//>, (107) 
which shows that a is determined by the residual between kinetic helicity of the small- 
scale velocity, (u? • u), and the normalized small-scale current density, (j ■ b)/p. 

The (j ■ b) term contributes to the nonlinear saturation of the dynamo. This is 
because the a effect produces magnetic helicity both at large and small scales such as 



to obey the magnetic helicity equation; see §5^ While this can lead to a dramatic 
reduction of a in periodic or closed domains |213[ I323j . the quenching effect may be less 
extreme in the astrophysically relevant case of open domains where magnetic helicity 
can be transported out of the domain by magnetic helicity fluxes |324[ 1325] . The theory 
of these fluxes |326j shows that there can be several contributions to the flux. One 
such contribution is along the contours of constant shear [3271 1328 j . but recent work has 
cast some doubt on whether such shear-driven magnetic helicity fluxes really exist 



Other contributions can come from advection [330J and diffusion [233J. For further 
aspects regarding nonlinear dynamo theory we refer to a review dedicated to recent 
developments; see Ref. [TT] . 

The presence of shear provides an additional induction effect that usually 
contributes to the dynamo. In order to estimate the relative importance of these effects, 



Astrophysical turbulence modeling 



59 



and to estimate whether a large-scale dynamo is excited, one needs to know the values 
of some relevant non-dimensional numbers that characterize the magnitude of a effect 
and shear, 

C a = a/ VT k m , C n = An/TfrfcL (108) 

where k m is an estimate for the relevant wavenumber of the dynamo that fits into 
the domain, and AQ is the absolute differential rotation. In the case of the Sun it is 
about 30% of the average angular velocity. Let us quantify the degree of helicity in the 
turbulence as e { = (u> ■ u)/kfU 2 ms , where kf = uj rras /u rras , and assume r] t rj, we find 

C a « e { ki/k m . (109) 

Thus, the efficiency of the a effect depends on how helical the turbulence is and on 
the amount of scale separation available. A so-called a 2 dynamo is possible when C a 
exceeds a critical value of the order of unity. 

Often C a is too small, and then the presence of shear can help to produce large- 
scale dynamo action. We assume that the shear is a significant fraction, q = AQ/Q, of 
the mean angular velocity Q which, in turn, is often expressed as the Coriolis number, 
Co = 2Q/u rras k { . We may then estimate 

C n ^lqCo(k { /k m ) 2 . (110) 

A large-scale dynamo of afl type is excited when the product C a Cn exceeds another 
critical value which is also of the order of unity. For a homogeneous dynamo, the critical 
value of C a Cn for plane wave solutions is 2. 

In conclusion, we see that the possibility of large-scale dynamo action depends 
critically on the scale separation ratio, kf/k m . It is therefore important that the domain 
is big enough to contain a significant number of turbulent eddies. Simulations have 
now confirmed the possibility of large-scale dynamo action in cases of forced turbulence, 
convective turbulence, and for turbulence driven in turn by magnetic fields through the 
magneto-rotational instability. 



7.3. Lambda effect 

An effect somewhat analogous to the a effect is the A effect. It parameterizes the 
dependence of the Reynolds stress on the mean angular velocity [33 U I332j as 

u^UTj = A ijk Q k +Afijki-^, (111) 

where Q = U^/ (r sin 9) is the local angular velocity (not the Qq used earlier in connection 
with the transformation to a rigidly rotating frame of reference). The second term in 



equation (111) is just the tensorial form of the turbulent viscosity; see equation (100). 
The first one exists already in the presence of uniform rotation. It is this term, balanced 
against the resulting turbulent viscosity term, that drives and maintains non-uniformity 
in the mean angular velocity [333, 1334[ 1335] . There are two important contributions to 
the A effect, a vertical and a horizontal one that quantify the r<p and 8(f) components 
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of the Reynolds stress, respectively. In particular, we have, in spherical coordinates, 
(r,M), 

/ Vsin#\ 

Atffcfifc = Fcos0 . (112) 

VVsintf Hcos6 / 

Here, V and if are functions of r and # that depend on the anisotropy of the turbulence. 
Using the first order smoothing approximation one finds |331[ I332j 

V = r(tf ) -u 2 ), (113) 
H = r(ul-uf). (114) 



For small turbulent Taylor numbers, Ta tur b = (2QR 2 /u t ) 2 , one finds for H = and 
that the Q contours are purely radial, while for V = and H ^ the Q contours 
are purely spoke-like. For V = z/ t (sin 2 # — 1) and H = u t sm 2 9 one finds disk-shaped fl 
contours. For V — Vt{\ sin 2 # — 1) and H = |z/ t sin 2 # one finds approximately spoke-like 
contours. However, those contours can change significantly with increasing values of 
Tat ur b, which leads to the development of cylindrical Q contours. This is explained by 
the Taylor- Proudman theorem, as will be explained below. 

The development of differential rotation is well established and is routinely seen in 
direct simulations of convective turbulence in rotating shells [551 1336, 337J . The existence 
of the A effect has been verified in local Cartesian simulations and the magnitude 
and spatial dependence has been determined |338j . Solutions of the equations for U 
have shown differential rotation roughly similar to what is found for the Sun using 
helioseismology. However, both DNS and solutions of the mean field equations show 
a tendency toward Q contours being constant along cylinders, which is not the case 
in the Sun. The cylindrical contours are the result of a feedback from the production 
of meridional circulation modifying the angular velocity contours. This leads to an 
approximate geostrophic balance, where 

U VU + p~ 1 Vp = 0. (115) 

In the barotropic case, when VT and Vs are parallel to each other, taking the curl 



of equation (115) yields V x (U ■ VE7) = 0. Assuming that the mean flow is purely 



toroidal, i.e. U = (0, 0, fir sin 6), we have 

on 2 

rsm6^— = 0. (116) 

oz 

So, if viscous and inertial terms are small, which is indeed the case for rapid rotation, 
then dVt /dz must be small, so f2 would be constant along cylinders |333j . This is also 
what is seen in mean-field models with A effect [339, 340J. 

It is generally believed that the main reason for Q not having cylindrical contours 
in the Sun is connected with the presence of the baroclinic term [? ]. The highest 
resolution simulations available today produce Q contours that are still too close to 
being constant along cylinders [S3 [551 H21[ 13411 13421 I343j 1344] . These simulations 
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Figure 25. Contours (left) and radial profiles (right) of differential rotation in a 
model of the Sun. Note the presence of radial deceleration near the surface layers 
(fractional radius above 0.9). Adapted from Ref. 



do not quite reach the solar surface, so they cannot show the near-surface shear layer 
where the rotation rate drops by more than 20nHz over the last 30 Mm just below the 
surface. Nevertheless, these simulations reproduce some important features of the Sun's 
differential rotation such as a more rapidly spinning equator. 

Mean-field simulations using the A effect show surprisingly good agreement with 
the helioseismologically inferred Q pattern |34U[ 1345] . and they are also beginning to 
reproduce the near-surface shear layer; see Figure |25| In these models it is indeed the 
baroclinic term that is responsible for causing the departure from cylindrical contours. 
This, in turn, is caused by an anisotropy of the turbulent heat conductivity which causes 
a slight enhancement in temperature and specific entropy at the poles. In the bulk of 
the convection zone the specific entropy is nearly constant while the temperature varies 
significantly in the vertical direction. It is therefore primarily the latitudinal specific 
entropy variation that determines the baroclinic term. This can be demonstrated by 
focusing on the contribution from the radial temperature and the latitudinal specific 
entropy gradients; see equations (49) and (53), i.e. 

QQ 2 „ i QT ds 

rsinfl— - « -0 ■ VT x Vs « — — — < 0. (117) 
oz r or 09 

The inequality shows that negative values of dfl jdz require that the pole is slightly 
warmer than the equator {ds/dd < 0). However, this effect is so weak that it cannot 
at present be observed. Allowing for these conditions in a simulation may require 
particular care in the treatment of the outer boundary conditions, or perhaps at the 
bottom of the convection zone in the tachocline. Given that the turbulent convective 
flux is proportional to — XijVjS, a negative ds/36 can be produced from a positive 
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enthalpy flux with a positive value of Xr6- This is indeed compatible with theory that 
predicts a rotational influence on the turbulence which makes Xij anisotropic |317[ |346j . 
One expects 

Xii = XtSij + Xne t , k nf } + X nM 0) nf\ (118) 

where we have used superscripts (0) interchangeably with subscripts to denote the 
rotation vector in a rotating frame of reference. In spherical polar coordinates we have 
O = (cos 9, — sin 9, 0), so Xre = ~ Xnn sin 9 cos 9 Qq- Simulations confirm that xnn is 
negative, but only when the scale of the mean field is comparable with that of the 
fluctuating velocity field |317j . which is somewhat unexpected. An alternative idea was 
advanced by Rempel |347j . who was able to reproduce solar- like fl contours by imposing 
a suitable latitudinal s variation at the bottom of the convection zone. 

In the discussion above we ignored in the last step the correlation between specific 
entropy and temperature fluctuations, i.e. a contribution from the term VX 7 x W 
where primes denote fluctuations. Such correlations, if of suitable sign, might provide 
yet another explanation for a non-zero value of dfl jdz. 

It is in principle also possible that the differential rotation could entirely be driven 
by the baroclinic term |348j 1349] . However, quantitative calculations showed that this 
effect on its own would be too small to |332[ [350J. 



7.4- Turbulent transport coefficients from simulations 

In the past few years significant progress has been made in determining tensor 
components such and r/ijk from local turbulence simulations. The 

recommended approach is what is referred to as the test-field method |351l 1352] . This 
method is not to be confused with the test-field model that was introduced by Kraichnan 
[353J as a closure approach. 



In the test-field method one solves numerically the evolution equation (92) for 
the fluctuations of the passive scalar concentration c, or a corresponding equation 
for fluctuations of the magnetic field b to obtain the magnetic transport coefficients. 
These equations are inhomogeneous in c or b and have terms of the form V • (uC) or 
V x (it x B). Here the mean fields C and B are now replaced by test fields. The best 
studied cases are for periodic boundary conditions and then the test fields are taken to 
be C° x = cos kx or C sx = sin kx, and similarly for the y and z directions. For each test 
field one evaluates the corresponding flux, UjC pq , and computes 

Kij = — (cos kxj J 13 * — smkxjT\ c ) /k, (119) 

for i,j = x,y,z. Here, angular brackets denote volume averages. Using this method it 
has now been possible to compute the dependence of the coefficients K t , kq, and kqq 



in an equation analogous to equation (118), where x h as been replaced by k. A similar 



equation can also be written down for the case where the anisotropy is caused by an 
applied magnetic field. 
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In the presence of a linear shear flow with Uij = const, it has proved advantageous 
to express «y in terms of the tensors Sy = \{U it j + Uj^) and Ay = \{Uij — Uj^). The 
corresponding representation of has been found to be of the form 

«ss(SS) ii + KAs(AS) ii . (120) 

There are indications that, in addition to K t , only the coefficients k$ and k$s are 
important, while ka and kas are either small or become small at larger Peclet number 

The test-field method also allows one to compute turbulent transport coefficients 
where the assumption of scale separation is not obeyed, or where the mean quantities 
vary on time scales comparable to the turnover time of the turbulence. In those cases 
we have to replace the multiplications in equations (100)-(102) by convolutions with 
integral kernels of the corresponding transport coefficients, e.g. 

U&(X, t) = - [ k[ C) (x,x', t, f) 9C ^2 t l d 3 x ' ( 12 1) 

J ox[ 
and likewise for the other equations. If the system is homogeneous and statistically 
stationary, the kernels depend only on the differences x — x' and t — t'. In such cases 
the convolution in real space becomes a multiplication in Fourier space. The test-field 
method yields directly the Fourier-transformed kernels if the test fields consist of sine 
and cosine functions |316j . and if they are made time-dependent [355J. By changing 
the wavenumber and/or the frequency of the test fields one can then obtain the full 
wavenumber and frequency dependence of the Fourier-transformed kernel functions that 
enables one to compute the kernels in real space via Fourier transformation. 

It turns out that, for a range of quite different physical circumstances, the k 
dependence can well be fitted to the form of a Lorentzian proportional to [l+(ak/ki) 2 ]~ 1 , 
where a is a fit parameter of the order of unity. The frequency dependence can be fitted 
to a function whose Fourier transform corresponds to a multiplicative contribution to 
the kernel of the form e~*/ T coswot, where loq and r -1 are fit parameters that are of the 
order of the inverse turnover time, u rms kf. 

This type of analysis has been carried out both for turbulent transport coefficients 
of both passive scalars and magnetic fields. In the magnetic case there is, in addition 
to the turbulent magnetic diffusivity tensor, also the a tensor |267] . By using test fields 
proportional to sine and cosine functions, one can compute a and rj t simultaneously. It 
turns out that in the kinematic regime a and rjt reach asymptotic values for ReM > 1 
where a ~ a and n t ~ rj t0 with 



a = |t(oj ■ u), r] W = |r(u 2 ), where r = (w rms A; f ) L . (122) 

Both show similar behavior as far as their wavenumber and frequency dependence is 
concerned |357j . 

Anisotropies of the flow yield anisotropics in the a and r\ tensors. In addition, since 
the magnetic field is an active vector, it backreacts on the flow and makes it anisotropic 
even if it was otherwise isotropic. It turns out that the a and 7] tensors are of the form 

ay = ctiSij + a 2 BiBj, rjij = rjxSij + ^BiBj, (123) 
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Figure 26. ReM-dependence of a and fj t . Both curves are normalized by a^. Adapted 
from Ref. [358] 



where we have assumed that only the unit vector of the mean field, B = B/\B\ is 
crucial. Simulations with the test-field method have revealed that in the quenched state 
«i and <xi have opposite sign [358J. The test-field method has been used to compute the 
dependence of r\ t and r] t on Re^ in the nonlinearly saturated state for B ~ £> eq , where 
B eq is the field strength where magnetic and kinetic energies are in equilibrium. It turns 
out that at large Rejv/ the effects of a and r/ t tend to balance each other. Furthermore, 
rjt does not show a sharp decline like Re^/, as would be the case in two dimensions, but, 
even though Re^ is already around 600, there remains a weak decrease of r] t , without 



any obvious indications that this trend might level off (Figure 26). 
When the a and rj tensors are multiplied by B, the result is [358 



CiijBj 



(124) 



where k m = fi Q J ■ B/B is an effective wavenumber of the mean field. This shows 
that the tensorial nature of a is unimportant in this context. However, this changes 
when considering passive vector equations that are equivalent to the induction equation, 
with a passive vector field B that is similar to the actual magnetic field, but it has 
no effect on the motions. Such a passive vector field can display dynamo action and 
can continue to grow even when the underlying velocity field corresponds to that of 
a nonlinearly saturated dynamo. This phenomenon was first observed for turbulent 
convective dynamos |359] and then confirmed for laminar dynamos generating a mean 



field that is an eigenvector of the matrix BiBj with vanishing eigenvalue. Thus, given 
that az/ai is negative, such fields remain unquenched for a velocity field or an a tensor 
that correspond to a saturated dynamo [360J. 
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8. Concluding remarks 

Over the past few decades hydrodynamic and magnetohydrodynamic simulations have 
become a frequently used tool in astrophysics research. This trend is surely going to 
continue. As an example of the importance of turbulence considerations we mention 
here the well-established field of stellar structure, which has recently been the subject 
of intense debate, because fits to three-dimensional time-dependent turbulent model 
atmospheres have led (mostly due to non-3D effects!) to a significantly lower estimate 
of the solar abundance of heavier elements. Although this issue is not yet settled, it 
is clear that the results from three-dimensional turbulence simulations will continue to 
provide valuable input to the debate. 

Even the radially symmetric (one-dimensional) models of stellar interiors are bound 
to be soon superseded or amended by higher-dimensional models. Clearly, the vast range 
of time and length scales between those of turbulent convection of stars and those of 
stellar structure and evolution necessitate a proper understanding of the collective or 
mean-field effects that are controlled by various correlators discussed in §[7| Obviously, 
we were only able to expose a small part of the many recent developments in this field. 
Quite frequently astrophysical turbulence involves magnetic fields, and often many more 
ingredients such as dust, chemicals, cosmic rays, and coupling to radiation. Instead of 
simply neglecting such additional features, one may attempt to incorporate them into 
stellar evolution models using a mean-field approach. The transport properties depend 
on rotation, shear, and magnetic field in ways that are reasonably well understood 
now. This is important for example in understanding the dependence of the Lithium 
abundance of young stars on their rotation rate |361j 362J. 

Astrophysical turbulence concerns usually extreme parameter regimes: large 
Reynolds and/or Mach numbers, very large or very small Prandtl numbers, as well 
as extreme density and temperature contrasts. This motivates thorough studies of 
turbulence in regimes that are not otherwise addressed. This can either provide broader 
support for certain turbulence theories, or it can more clearly highlight problems that 
would be otherwise overlooked. In this sense astrophysical turbulence research is not 
just the application of regular turbulence theory, but it can also provide complementary 
insights of broader relevance also for other research fields. 

One of the aspects where astrophysical turbulence encounters an as yet unsettled 
issue is the question how compressibility really enters the theory. We have seen some 
ambiguity in the proper definition of the kinetic energy where, empirically, the spectrum 
of p x l z u appears to be closest to the case of incompressible turbulence. There are several 
related issues in the context of mean-field theory. For example, the equation for the 



magnetic a effect in equation (107) contains a p factor, but since this equation was 
derived for the compressible case, it is not clear whether p should enter inside or outside 
the average of j ■ b when p is non-uniform or strongly fluctuating. Another occurrence of 
compressibility effects could be in the relation between the enthalpy flux and the specific 
entropy gradient. Finally, let us also mention here the issue of the baroclinic term, which 
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can be important in the production of vorticity and shaping the form of the differential 
rotation contours in the Sun. There could potentially be systematic corrections resulting 
from the fluctuations of specific entropy and temperature. This and other effects might 
be responsible for causing a departure from cylindrical Taylor-Proudman contours of 
Q(r, 6) in the Sun. 

There are several other quadratic correlation functions that need to be modelled 
more accurately. One is the current helicity flux involving terms of the form E x J, for 
example. Other examples include Reynolds and Maxwell stresses and their dependence 
not only on the mean velocity, as discussed above, but also on the magnetic field. This 
quadratic nonlinearity means that the standard test-field method cannot be used, but 
possibly some kind of modification of it might work. 
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